Bayesian Learning-Based Clustered-Sparse Channel Estimation for Time-Varying Underwater Acoustic OFDM Communication ^†

Abstract

Orthogonal frequency division multiplexing (OFDM) has been widely adopted in underwater acoustic (UWA) communication due to its good anti-multipath performance and high spectral efficiency. For UWA-OFDM systems, channel state information (CSI) is essential for channel equalization and adaptive transmission, which can significantly affect the reliability and throughput. However, the time-varying UWA channel is difficult to estimate because of excessive delay spread and complex noise distribution. To this end, a novel Bayesian learning-based channel estimation architecture is proposed for UWA-OFDM systems. A clustered-sparse channel distribution model and a noise-resistant channel measurement model are constructed, and the model hyperparameters are iteratively optimized to obtain accurate Bayesian channel estimation. Accordingly, to obtain the clustered-sparse distribution, a partition-based clustered-sparse Bayesian learning (PB-CSBL) algorithm was designed. In order to lessen the effect of strong colored noise, a noise-corrected clustered-sparse channel estimation (NC-CSCE) algorithm was proposed to improve the estimation accuracy. Numerical simulations and lake trials are conducted to verify the effectiveness of the algorithms. Results show that the proposed algorithms achieve higher channel estimation accuracy and lower bit error rate (BER).

1. Introduction

Underwater acoustic (UWA) sensor networks can be applied to a variety of underwater activities, including environment monitoring, resource exploration and equipment navigation. Reliable high-rate communication is significant for massive data transmission between sensor nodes [1]. However, the UWA channel is widely known as one of the most difficult communication media, and high-rate communication with high reliability is quite challenging. Since the sound attenuation rapidly increases with frequency, the bandwidth available for communication is extremely limited. Due to the reflection and refraction of sound in the ocean, the multipath spread in a medium-range shallow water channel can extend over tens or even hundreds of milliseconds [2], resulting in severe frequency-selective signal distortion. Due to the dynamic environment and moving transceivers in ocean, the UWA channel is time-varying. In particular, the normalized frequency offset factor induced by Doppler can be on the order of ${10}^{-3}$ compared with ${10}^{-7}$ for high-speed mobile radio channels [2].

Due to high spectral efficiency and good anti-multipath performance, orthogonal frequency division multiplexing (OFDM) and some OFDM-based modulation techniques are widely adopted for UWA communications [3,4,5,6]. In UWA-OFDM systems, accurate channel state information (CSI) is vital for channel equalization and adaptive transmission [7,8], which can significantly affect the reliability and throughput. However, it is difficult to estimate the time-varying channel impulse responses (CIRs) with excessive delay spread. To overcome this problem, some structural information such as sparsity can be exploited to improve channel estimation performance. The paths of the UWA channel tend to be clustered-sparse, which carry comprehensive structural information [9,10,11,12,13]. As shown in Figure 1, there are many weak paths centering around the eigen-paths associated with the medium refractions and surface/bottom bounces, which form the clustered-sparse structure. As different clusters are usually derived from different transmission routes, the channel coherence time usually varies across clusters. Thus, the path variation for different clusters can reflect the channel variation characteristics [14], and can be used for predicting the time-varying UWA channel [15]. However, the fading process of the UWA channel may randomly vary [2], which makes it difficult to accurately obtain the path variation for different clusters between adjacent OFDM blocks. To obtain the clustered-sparse structure of the UWA channel, we proposed a clustered-sparse Bayesian channel estimation method in [16], which is the preliminary version of this work (This work extends the conference paper [16]) by (i) adding a Bayesian learning-based channel estimation architecture; (ii) adding a noise-corrected clustered-sparse channel estimation algorithm to lessen the effect of strong colored noise; (iii) giving the detailed derivation and performance analysis of the proposed algorithms; and (iv) presenting more complete performance results related to modulation order, model hyperparameters and noise distribution).

In addition to the time-varying CIRs with excessive delay spread, the complex distribution of noise in the UWA channel also makes the channel estimation difficult. Generally, there are two kinds of channel noise, including the ambient noise and site-specific noise. The former type is presented in the background of the quiet deep sea, and exhibits a colored characteristic when approximated by Gaussian distribution, while for the site-specific noise, significant non-Gaussian components should be carefully considered [2]. Some experiments have confirmed that the noise is non-uniform across different locations, times and frequencies [17]. As a popular noise model in most UWA-OFDM systems, additive white Gaussian noise (AWGN) is straightforward for modeling and solving the noise distribution of the UWA channel, but it deviates from the actual underwater noise. Obviously, the deviation of the noise model means the estimation error of the noise distribution, which will lead to the degradation of channel estimation and signal detection performance.

To tackle the challenge of estimating the time-varying UWA channel with excessive delay spread and complex noise distribution, we adopted Bayesian learning to obtain the multipath and noise distribution characteristics of the UWA channel, based on which Bayesian channel estimation is derived to improve the performance of pilot-based channel estimation. The main contributions of this work are summarized as follows:

• To estimate the time-varying multipath channel with colored noise, we propose a novel Bayesian learning-based channel estimation architecture for UWA-OFDM systems. Specifically, a clustered-sparse channel distribution model is constructed to characterize the delay power spectrum and temporal correlation of each cluster in the multipath channel, and a noise-resistant channel measurement model is constructed to reduce the noise disturbance. By learning the model hyperparameters, the Bayesian channel estimation based on the two models can be iteratively optimized.
• To obtain the clustered-sparse distribution, we propose a partition-based clustered-sparse Bayesian learning (PB-CSBL) algorithm. Through the cluster partition, different clusters can learn different channel correlation coefficients, and thus the inter-cluster interference of the multipath channel can be suppressed.
• To lessen the effect of strong colored noise, we propose a noise-corrected clustered-sparse channel estimation (NC-CSCE) algorithm. Based on the iterative symbol decision and noise correction, the more accurate hyperparameters of the models can be obtained, which can improve the accuracy of the Bayesian channel estimation.

The rest of this paper is organized as follows. Section 2 discusses the related works. In Section 3, the system architecture and channel models are constructed. Then, the Bayesian learning-based clustered-sparse channel estimation method is designed in Section 4. In Section 5, the noise-corrected clustered-sparse channel estimation method is presented. The evaluation and result analysis are shown in Section 6. In Section 7, conclusions are given.

2. Related Works

The UWA channel is severely band-limited, and exhibits time-varying excessive delay spread. To this end, the sparse structure of the UWA channel can be exploited to improve the estimation performance, as most channel energy of the UWA channel is concentrated on a few paths [18,19]. The paths of the UWA channel also tend to be clustered-sparse [20], which has been frequently observed based on field experiments at different locations [9,10,11,12,13]. A detailed investigation on the variation of the significant paths in each cluster is presented in [21]. Based on the observation of channel sparsity, compressed sensing (CS)-based sparse channel estimation methods [22,23,24] have been studied, such as the greedy matching pursuit (MP) and orthogonal matching pursuit (OMP), which were adopted in [18,25], aiming to estimate the UWA channel impulse response (CIR) and the channel delay-Doppler-spread function, respectively. In order to reduce the convergence error of the OMP algorithm, the compressed sampling matching pursuit (CoSaMP) algorithm is proposed in [26]. However, CoSaMP-based channel estimators require knowledge of the channel sparsity level, which is difficult to obtain accurately in advance in practical applications. The joint sparsity among adjacent OFDM blocks was exploited by the simultaneous OMP (SOMP) [27] for higher estimation accuracy, but only stable path delays were considered and the temporal correlation of path gains was ignored. Typical convex optimization algorithms including basis pursuit (BP) and basis pursuit denoising (BPDN) tackle the NP-hard $l_0$ norm problem by transforming the $l_0$ norm to a convex $l_1$ norm. The BP algorithm was applied to estimate the path delays and Doppler scale factors of the UWA channel in [25], and the BPDN algorithm was proposed to estimate the slowly time-varying UWA channels in [28]. However, the global minimum of the cost function in BP or BPDN is not necessarily consistent with the sparsest solution. An empirical channel variation model, in which different clusters vary independently, was proposed in [14], and demonstrated that adaptation with multiple clusters outperforms that with one cluster for UWA channel estimation.

In recent years, Bayesian learning (BL) has been used to improve the performance of channel estimation. Compared with the CS-based methods that provide point estimation, Bayesian techniques can give the interval estimation of channel response and have a desirable property of preventing structural errors [29]. By using BL-based methods, the high accuracy of channel estimation was achieved in OFDM systems [30,31,32,33]. The temporal multiple sparse Bayesian learning (TMSBL)-based estimation method was further proposed in [32] to estimate the sparse channels by taking advantage of the channel coherence between consecutive OFDM blocks, and it can achieve better performance in strong correlated channels and maintains robustness in weak temporal correlated channels. Due to the high computational complexity of BL, some accelerated algorithms were proposed, such as the approximate message passing (AMP)-based method in [34]. More recently, some researchers utilized deep learning (DL)-based channel estimators for OFDM systems [35,36,37], which showed better channel estimation or equalization performance for nonlinear channel distortion. However, how to go about improving the generalization performance of pre-training DL models and reduce the complexity of training and inferring remain challenging.

3. System Architecture and Channel Models

In this section, the system architecture of Bayesian learning-based channel estimation is presented first. Then, the noise-resistant channel measurement model and the clustered-sparse channel distribution model are constructed.

3.1. Bayesian Learning-Based Channel Estimation Architecture

A typical UWA-OFDM communication scenario was shown at the top of Figure 2. The receiver can obtain the transmitted OFDM signal through dynamic multipaths, and the received signal may be affected by random colored noise. At the bottom of Figure 2, the system architecture of Bayesian learning-based channel estimation is given. To characterize the delay power spectrum and temporal correlation of each cluster in the UWA channel, the clustered-sparse channel distribution model is constructed and the model hyperparameters can be optimized by the iterative partition-based cluster evolution. Demodulated signal and pilot symbols were utilized to measure the instantaneous channel response, and the iterative noise measurement and data detection can improve the anti-noise performance of the channel measurement model. Based on the channel distribution model and the channel measurement model, the Bayesian channel estimation can be obtained.

Important notations used in following sections are given by:

• Uppercase and lowercase bold symbols are reserved for matrices and vectors, respectively. Particularly, ${\mathbf{I}}_{\mathbf{M}}$ denotes the identity matrix with size $M\times M$. When the dimension is evident from the context, for simplicity, we only used $\mathbf{I}$;
• ${\mathbf{A}}^T$ and ${\mathbf{A}}^H$ represent the transpose and conjugate transpose of $\mathbf{A}$, respectively;
• $\mathbf{A}$⊗$\mathbf{B}$ denotes the Kronecker product of the two matrices $\mathbf{A}$ and $\mathbf{B}$;
• $vec\left(\mathbf{A}\right)$ denotes the vectorization of $\mathbf{A}$ formed by stacking its columns into a single column vector;
• $Tr\left(\mathbf{A}\right)$ denotes the trace of $\mathbf{A}$;
• $Toeplitz(a_1,a_2,\dots ,a_N)$ denotes a Toeplitz matrix taking $a_1,\dots ,a_N$ as first row;
• $diag\left(\mathbf{a}\right)$ or $diag(a_1,a_2,\dots ,a_N)$ denotes a diagonal matrix with principal diagonal elements being $a_1,a_2,\dots ,a_N$ in turn; $diag({\mathbf{A}}_1,{\mathbf{A}}_2,\dots ,{\mathbf{A}}_N)$ denotes a block diagonal matrix with principal diagonal blocks being the square matrices ${\mathbf{A}}_1,\dots ,{\mathbf{A}}_N$ in turn;
• If $\mathbf{A}$ is a square matrix, $Prdiag\left(\mathbf{A}\right)$ denotes a diagonal matrix with principal diagonal elements being the principal diagonal elements of $\mathbf{A}$ in turn.
• If some terms in a cost function do not contribute to the subsequent optimization of the parameters, ∝ is used to indicate that these terms have been dropped.

3.2. Noise-Resistant Channel Measurement Model

Considering the OFDM system with cyclic prefix (CP) and K subcarriers, the non-overlapping sets of $K_d$ data subcarriers ${\mathbf{x}}_{\mathbf{d}}$, $K_p$ pilot subcarriers ${\mathbf{x}}_{\mathbf{p}}$ and $K_n$ null subcarriers ${\mathbf{x}}_{\mathbf{n}}$ satisfy $K=K_d+K_p+K_n$. Let T denotes the OFDM symbol period without CP, $T_{cp}$ denotes the length of CP and $T_s=T+T_{cp}$ is the whole OFDM symbol period. The kth subcarrier is at the frequency:

$$f_k=f_c+k/T,k=-K/2,\dots ,K/2-1,$$

(1)

where $f_c$ is the carrier frequency. Assuming that there are M OFDM blocks in one frame, a passband waveform of the mth OFDM block at time $t^{\prime }\in [-T_{cp}+m{T}_s,T+m{T}_s]$ is then given by

$${\tilde{x}}^{\left(m\right)}\left(t\right)=2Re\left\{\frac{1}{\sqrt{K}}\sum _{k=-K/2}^{K/2-1}{x}_k^{\left(m\right)}{e}^{j2\pi f_{k}t}q\left(t\right)\right\},$$

(2)

where $m=1,2,\dots ,M$, $t=t\prime -m{T}_s$, $x_k^{\left(m\right)}$ is defined as a transmitted symbol on the kth subcarrier at the mth OFDM block and q(t) is the pulse shaping filter:

$$q\left(t\right)=\left\{\begin{array}{cc}\hfill 1& ,t\in [-T_{cp},T],\hfill \\ \hfill 0& ,otherwise.\hfill \end{array}\right.$$

(3)

The UWA channel is a typical sparse multipath channel. For horizontal shallow water multi-path transmission with a range much greater than depth, each path has almost the same Doppler factor [38]. Thus, we assume that all the paths have the equal Doppler scale factor a, and the Doppler effect can be almost eliminated through effective Doppler estimation and compensation. Furthermore, the path delay remains stable across several consecutive OFDM blocks, and the path gains and Doppler scale factors are constant during one OFDM block, but vary from block to block. The time-varying channel response during the mth OFDM block can be expressed as

$$h^{\left(m\right)}(\tau ,t)=\sum _{l=0}^{L-1}{A}_l^{\left(m\right)}\delta \left[\tau -\left({\tau }_l-a^{\left(m\right)}t\right)\right],$$

(4)

where L is the path number of the multipath channel, which may be distributed in more than one cluster. $A_l$ and ${\tau }_l$ denote the gain and delay of the lth path.

Through the time-varying channel, the received passband signal during the mth OFDM block can be expressed as

$${\tilde{y}}^{\left(m\right)}\left(t\right)=\sum _{l=0}^{L-1}{A}_l^{\left(m\right)}{\tilde{x}}^{\left(m\right)}\left[\left(1+a^{\left(m\right)}\right)t-{\tau }_l\right]+{\tilde{w}}^{\left(m\right)}\left(t\right),$$

(5)

where $\tilde{w}\left(t\right)$ is the additive noise. After synchronization, a popular two-step approach [39] can be adopted to mitigate the Doppler effect, which takes a coarse Doppler estimation and resampling as the first step and then performs the fine Doppler shift compensation using the null subcarriers. Moreover, CP can also be used for Doppler shift estimation [40]. Performing CP-OFDM demodulation, the output $K\times 1$ vector of the mth OFDM block can be expressed as

$${\mathbf{y}}^{\left(m\right)}={\mathbf{X}}^{\left(m\right)}\mathbf{F}{\mathbf{h}}^{\left(m\right)}+{\mathbf{w}}^{\left(m\right)},$$

(6)

where ${\mathbf{X}}^{\left(m\right)}$ is an $K\times K$ diagonal matrix whose diagonal entries are the K transmitted symbols, $\mathbf{F}$ is an $K\times L$ discrete Fourier transform (DFT) matrix and ${\mathbf{w}}^{\left(m\right)}$ is an $K\times 1$ additive noise vector. ${\mathbf{h}}^{\left(m\right)}$ is the discretized CIR vector as ${\mathbf{h}}^{\left(m\right)}=[h_1^{\left(m\right)},h_2^{\left(m\right)},\dots ,h_L^{\left(m\right)}]$.

For M consecutive OFDM blocks, the transmission model of all $K\times M$ symbols can be expressed as

$$\mathbf{Y}=\left[{\mathsf{\Phi }}^{\left(1\right)}{\mathbf{h}}^{\left(1\right)},{\mathsf{\Phi }}^{\left(2\right)}{\mathbf{h}}^{\left(2\right)},\dots ,{\mathsf{\Phi }}^{\left(M\right)}{\mathbf{h}}^{\left(M\right)}\right]+\mathbf{W},$$

(7)

where $\mathbf{Y}=[{\mathbf{y}}^{\left(1\right)},{\mathbf{y}}^{\left(2\right)},\dots ,{\mathbf{y}}^{\left(M\right)}]$, $\mathbf{W}=[{\mathbf{w}}^{\left(1\right)},{\mathbf{w}}^{\left(2\right)},\dots ,{\mathbf{w}}^{\left(M\right)}]$, and $\mathbf{H}=[{\mathbf{h}}^{\left(1\right)},{\mathbf{h}}^{\left(2\right)},\dots ,{\mathbf{h}}^{\left(M\right)}]$. ${\mathsf{\Phi }}^{\left(\mathbf{m}\right)}$ is defined as ${\mathsf{\Phi }}^{\left(\mathbf{m}\right)}\triangleq {\mathbf{X}}^{\left(m\right)}\mathbf{F}$. Considering the comb-type pilot arrangement, the transmission model of $K_p\times M$ pilot symbols can be expressed as

$${\mathbf{Y}}_{\mathbf{P}}={\mathsf{\Phi }}_{\mathbf{P}}\mathbf{H}+{\mathbf{W}}_{\mathbf{P}},$$

(8)

where ${\mathbf{Y}}_{\mathbf{P}}=[{\mathbf{y}}_{\mathbf{p}}^{\left(1\right)},{\mathbf{y}}_{\mathbf{p}}^{\left(2\right)},\dots ,{\mathbf{y}}_{\mathbf{p}}^{\left(M\right)}]$ and ${\mathbf{W}}_{\mathbf{P}}=[{\mathbf{w}}_{\mathbf{p}}^{\left(1\right)},{\mathbf{w}}_{\mathbf{p}}^{\left(2\right)},\dots ,{\mathbf{w}}_{\mathbf{p}}^{\left(M\right)}]$ are the submatrixs of $\mathbf{Y}$ and $\mathbf{W}$ at the location of the pilot subcarriers, respectively. ${\mathsf{\Phi }}_{\mathbf{P}}\triangleq {\mathbf{X}}_{\mathbf{P}}{\mathbf{F}}_{\mathbf{P}}$ is the known dictionary matrix where ${\mathbf{X}}_{\mathbf{P}}$ is the $K_p\times K_p$ diagonal matrix with the known pilots along its diagonal and ${\mathbf{F}}_{\mathbf{P}}$ is the $K_p\times L$ DFT matrix.

By vectorizing $\mathbf{H}$ as $\mathbf{h}=vec\left({\mathbf{H}}^T\right)$, a pilot-based channel measurement model can be constructed according to (8), as

$${\mathbf{y}}_{\mathbf{p}}={\mathbf{U}}_{\mathbf{P}}\mathbf{h}+{\mathbf{w}}_{\mathbf{p}},$$

(9)

where ${\mathbf{y}}_{\mathbf{p}}=vec\left({\mathbf{Y}}_{\mathbf{P}}^T\right)$, ${\mathbf{w}}_{\mathbf{p}}=vec\left({\mathbf{W}}_{\mathbf{P}}^T\right)$ and ${\mathbf{U}}_{\mathbf{P}}={\mathsf{\Phi }}_{\mathbf{P}}\otimes {\mathbf{I}}_M$. AWGN can be adopted to approximate the additive noise ${\mathbf{w}}_{\mathbf{p}}$. However, the attenuation and noise of the UWA channel may vary over the signal bandwidth. When a small number of pilot symbols are used to estimate the UWA channel with strong colored noise, the accuracy of the pilot-based channel measurement will decrease severely. To improve the performance of channel measurement, a noise-resistant channel measurement model is constructed according to (7) and it can be expressed as

$$\mathbf{y}=\mathbf{U}\mathbf{h}+\mathbf{w},$$

(10)

where $\mathbf{y}=vec\left({\mathbf{Y}}^T\right)$, $\mathbf{w}=vec\left({\mathbf{W}}^T\right)$ and $\mathbf{U}$ is the $KM\times LM$ matrix as

$$\mathbf{U}=\left[\begin{array}{ccc}diag\{{\mathsf{\Phi }}^{\left(1\right)}[1,1],\dots ,{\mathsf{\Phi }}^{\left(M\right)}[1,1]\}& \dots & diag\{{\mathsf{\Phi }}^{\left(1\right)}[1,L],\dots ,{\mathsf{\Phi }}^{\left(M\right)}[1,L]\}\\ ⋮& \ddots & ⋮\\ diag\{{\mathsf{\Phi }}^{\left(1\right)}[K,1],\dots ,{\mathsf{\Phi }}^{\left(M\right)}[K,1]\}& \dots & diag\{{\mathsf{\Phi }}^{\left(1\right)}[K,L],\dots ,{\mathsf{\Phi }}^{\left(M\right)}[K,L]\}\end{array}\right],$$

(11)

where ${\mathsf{\Phi }}^{\left(m\right)}[i,j]$ denotes the entry in the ith row and the jth column of ${\mathsf{\Phi }}^{\left(m\right)}$. All subcarriers are used for channel measurement in (10), as the unknown data symbols can be obtained by data detection or decoding. Assuming that the power distribution over all subcarriers is non-uniform, $\lambda =[{\lambda }_{-K/2},{\lambda }_{-K/2+1},\dots ,{\lambda }_{K/2-1}]$, $\mathbf{w}$ in (10) satisfies:

$$\mathbf{w}\sim \mathcal{CN}(\mathbf{0},\mathsf{\Lambda }),$$

(12)

where $\mathsf{\Lambda }=diag\left(\lambda \right)\otimes {\mathbf{I}}_M$.

3.3. Clustered-Sparse Channel Distribution Model

According to the measurement equation in (9) or (10), the channel vector $\mathbf{h}$ can be estimated. However, the UWA channel usually has an excessive multipath delay spread. With so many undetermined channel taps, the prior clustered-sparse structure of the UWA channel can be exploited to obtain the sparse solution of the channel vector. First, we assume all that the L paths (i.e., L rows) in $\mathbf{H}$ are mutually independent, which can form C clusters. There are $L_d$ paths in the dth cluster and these clusters do not overlap, each of which occupies dense path delays. With the same time correlation characteristics, the dth cluster satisfies a first-order auto-regressive (AR) model, given by

$${\mathbf{H}}_d[i,j+1]={\beta }_d{\mathbf{H}}_d[i,j]+\sqrt{1-|{\beta }_d{|}^2}{\mathbf{V}}_d[i,j],$$

(13)

where ${\beta }_d$ is the temporal correlation coefficient of the dth cluster. ${\mathbf{H}}_d$ is a $L_d\times M$ submatrix of $\mathbf{H}$, $i=1,2,\dots ,L_d$ and $j=1,2,\dots ,M-1$ are the row and column indices of ${\mathbf{H}}_d$, respectively. ${\mathbf{V}}_d$ is the model noise matrix with complex Gaussian distribution and assumed that $vec\left({\mathbf{V}}_d^T\right)\sim \mathcal{CN}(0,{\mathsf{\Gamma }}_d\otimes {\mathbf{I}}_M)$, where:

$${\mathsf{\Gamma }}_d=diag\left({\gamma }_{d,1},{\gamma }_{d,2},\dots ,{\gamma }_{d,L_d}\right)$$

(14)

is a positive semi-definite diagonal real-valued matrix. Although a higher-order AR model may express better model approximation performance, it would lead to higher complexity. More importantly, a higher-order AR model means a higher overfitting risk, especially when there are only a few OFDM blocks in one frame. Hence, a first-order AR process was adopted to model the temporal correlation characteristic for each cluster of CIRs.

By vectorizing ${\mathbf{H}}_{\mathbf{d}}$ as ${\mathbf{h}}_{\mathbf{d}}=vec\left({\mathbf{H}}_d^T\right)$, the parameteric Bayesian prior of ${\mathbf{h}}_d$ can be written by

$$p({\mathbf{h}}_d;{\beta }_d,{\gamma }_d)=\mathcal{CN}(\mathbf{0},{\mathsf{\Gamma }}_d\otimes {\mathbf{B}}_d),$$

(15)

where ${\gamma }_d=\left[{\gamma }_{d,1},{\gamma }_{d,2},\dots ,{\gamma }_{d,L_d}\right]$ controls the sparsity of the $L_d$ paths in the dth cluster. When $Tr\left({\mathsf{\Gamma }}_d\right)\to 0$, the associated ${\mathbf{h}}_d\to \mathbf{0}$. In other words, ${\mathsf{\Gamma }}_d$ reflects the probability density of the dth cluster. ${\mathbf{B}}_d$ is a Toeplitz matrix with the following form:

$${\mathbf{B}}_d=Toeplitz\left(1,{\beta }_d,\dots ,{\beta }_d^{M-1}\right),$$

(16)

and it reflects the variation characteristics of the dth cluster. ${\mathsf{\Gamma }}_d$ and ${\mathbf{B}}_d$ jointly determine the cluster distribution of the dth cluster. Accordingly, the prior probability density function of $\mathbf{h}$ can be written as

$$p(\mathbf{h};{\gamma }_1,{\gamma }_2,\dots ,{\gamma }_C,{\beta }_1,{\beta }_2,\dots ,{\beta }_C)=\mathcal{CN}(\mathbf{0},{\sum }_0),$$

(17)

where ${\sum }_0=diag\left({\mathsf{\Gamma }}_1\otimes {\mathbf{B}}_1,{\mathsf{\Gamma }}_2\otimes {\mathbf{B}}_2,\dots ,{\mathsf{\Gamma }}_C\otimes {\mathbf{B}}_C\right)$.

4. Bayesian Learning-Based Clustered-Sparse Channel Estimation

In this section, the parametric form of Bayesian channel estimation is given first, and followed by the cluster distribution learning. Then, the partition-based clustered-sparse Bayesian learning algorithm is presented to obtain the clustered-sparse distribution and Bayesian channel estimate. At last, the complexity and performance analysis of this algorithm is given.

4.1. Bayesian Channel Estimation

To simplify the implementation complexity of Bayesian channel estimation, only known pilot symbols are used for channel measurement. As pilot symbols are sparsely distributed, it is difficult to estimate the non-uniform power spectrum at the location of pilot subcarriers. Therefore, the additive noise in (9) is approximated as white Gaussian noise with noise variance $\lambda$ in this section. According to (9) and (17), using the Bayes rule, the posterior density of $\mathbf{h}$ can be obtained as

$$p\left(\mathbf{h}\right|{\mathbf{y}}_{\mathbf{p}};\lambda ,{\gamma }_1,{\gamma }_2\dots ,{\gamma }_C,{\beta }_1,{\beta }_2,\dots {\beta }_C)=\mathcal{CN}({\mu }_{\mathbf{h}},{\sum }_{\mathbf{h}}),$$

(18)

which is the complex Gaussian distribution with the posterior mean and covariance matrix:

$${\mu }_{\mathbf{h}}={\lambda }^{-1}{\sum }_{\mathbf{h}}{\mathbf{U}}_{\mathbf{P}}^H{\mathbf{y}}_{\mathbf{p}}={\sum }_0{\mathbf{U}}_{\mathbf{P}}^H{\left(\lambda \mathbf{I}+{\mathbf{U}}_{\mathbf{P}}{\sum }_0{\mathbf{U}}_{\mathbf{P}}^H\right)}^{-1}{\mathbf{y}}_{\mathbf{p}},$$

(19)

$${\sum }_{\mathbf{h}}={\left({\sum }_0^{-1}+{\lambda }^{-1}{\mathbf{U}}_{\mathbf{P}}^H{\mathbf{U}}_{\mathbf{P}}\right)}^{-1}={\sum }_0-{\sum }_0{\mathbf{U}}_{\mathbf{P}}^H{\left(\lambda \mathbf{I}+{\mathbf{U}}_{\mathbf{P}}{\sum }_0{\mathbf{U}}_{\mathbf{P}}^H\right)}^{-1}{\mathbf{U}}_{\mathbf{P}}{\sum }_0.$$

(20)

Obviously, the maximum posterior (MAP) estimates ${\widehat{\mathbf{h}}}_{MAP}\triangleq {\mu }_{\mathbf{h}}$ and the covariance matrix ${\sum }_{\mathbf{h}}$ are related to the hyperparameters $\left[\lambda ,{\gamma }_1,{\gamma }_2,\dots ,{\gamma }_C,{\beta }_1,{\beta }_2,\dots ,{\beta }_C\right]$. Using $K_n\times M$ null subcarriers ${\mathbf{Y}}_{\mathbf{N}}=[{\mathbf{y}}_n^{\left(1\right)},\dots ,{\mathbf{y}}_n^{\left(M\right)}]$, $\lambda$ can be estimated as

$$\lambda =\frac{1}{K}E[\parallel {\mathbf{y}}^{\left(m\right)}{\parallel }_2^2]\approx \frac{1}{M{K}_n}\sum _{j=1}^M\sum _{i=1}^{K_n}{\left|{\mathbf{Y}}_{\mathbf{N}}[i,j]\right|}^2,$$

(21)

where i and j are the row and column indexes of ${\mathbf{Y}}_{\mathbf{N}}$, respectively. The remaining hyperparameters $\left[{\gamma }_1,{\gamma }_2,\dots ,{\gamma }_C,{\beta }_1,{\beta }_2,\dots ,{\beta }_C\right]$ are related to the cluster distribution information.

4.2. Cluster Distribution Learning

The hyperparameters $\mathbf{\theta }\triangleq \left[{\gamma }_1,{\gamma }_2,\dots ,{\gamma }_C,{\beta }_1,{\beta }_2,\dots ,{\beta }_C\right]$ can be estimated by maximizing the marginal likelihood function $p({\mathbf{y}}_{\mathbf{p}};\mathbf{\theta })$. This is equivalent to minimizing $-logp({\mathbf{y}}_{\mathbf{p}};\mathbf{\theta })$, giving the effective cost function:

$$\mathcal{L}\left(\mathbf{\theta }\right)={\mathbf{y}}_{\mathbf{p}}^H{\sum }_{\mathbf{y}}{\mathbf{p}}^{-1}\left(\mathbf{\theta }\right){\mathbf{y}}_{\mathbf{p}}+log\left|{\sum }_{\mathbf{y}}\left(\mathbf{\theta }\right)\right|,$$

(22)

where ${\sum }_{\mathbf{y}}\mathbf{p}\left(\mathbf{\theta }\right)=\lambda \mathbf{I}+{\mathbf{U}}_{\mathbf{P}}{\sum }_0\left(\mathbf{\theta }\right){\mathbf{U}}_{\mathbf{P}}^H$. The above problem cannot be solved in closed form, and the expectation maximization (EM) method [41] can be employed to solve it iteratively. For unknown values of hyperparameters governing the prior density, $\mathbf{h}$ is considered as the nuisance variable and $\mathbf{\theta }$ is estimated by maximizing:

$$Q\left(\mathbf{\theta }\right|{\mathbf{\theta }}^{\left(old\right)})=E_{\mathbf{h}|{\mathbf{y}}_{\mathbf{p}};{\mathbf{\theta }}^{\left(old\right)}}\left[logp({\mathbf{y}}_{\mathbf{p}},\mathbf{h};\mathbf{\theta })\right]=E_{\mathbf{h}|{\mathbf{y}}_{\mathbf{p}};{\mathbf{\theta }}^{\left(old\right)}}\left[logp\left({\mathbf{y}}_{\mathbf{p}}\right|\mathbf{h};\lambda )\right]+E_{\mathbf{h}|{\mathbf{y}}_{\mathbf{p}};{\mathbf{\theta }}^{\left(old\right)}}\left[logp(\mathbf{h};\mathbf{\theta })\right],$$

(23)

where ${\mathbf{\theta }}^{\left(old\right)}$ denotes the estimated hyperparameters in the previous iteration. Ignoring the first term unrelated to $\mathbf{\theta }$, the Q Function (23) can be simplified to:

$$\begin{array}{cc}\hfill Q({\beta }_1,{\beta }_2,\dots ,{\beta }_C,& {\gamma }_1,{\gamma }_2,\dots ,{\gamma }_C)=E_{\mathbf{h}|{\mathbf{y}}_{\mathbf{p}};{\mathbf{\theta }}^{\left(old\right)}}\left\{logp\left[\mathbf{h};{\mathbf{B}}_1,{\mathbf{B}}_2,\dots ,{\mathbf{B}}_C,{\mathsf{\Gamma }}_1,{\mathsf{\Gamma }}_2,\dots ,{\mathsf{\Gamma }}_C\right]\right\}\hfill \\ \hfill \propto & -\frac{1}{2}\sum _{d=1}^C\left[log\left(\right|{\mathsf{\Gamma }}_d{|}^M|{\mathbf{B}}_d{|}^{L_d})+{\mathbf{h}}_d^H({\mathsf{\Gamma }}_d^{-1}\otimes {\mathbf{B}}_d^{-1}){\mathbf{h}}_d\right]\hfill \\ \hfill \propto & -\frac{1}{2}\sum _{d=1}^C\left\{Mlog\left(\right|{\mathsf{\Gamma }}_d\left|\right)+L_{d}log\left(\right|{\mathbf{B}}_d\left|\right)+Tr\left[\left({\mathsf{\Gamma }}_d^{-1}\otimes {\mathbf{B}}_d^{-1}\right)\left({\sum }_{\mathbf{h}}^d+{\mu }_{\mathbf{h}}^d{\left({\mu }_{\mathbf{h}}^d\right)}^H\right)\right]\right\}\hfill \end{array},$$

(24)

where ${\mathbf{h}}_d\in {\mathbb{C}}^{M{L}_d\times 1}$, ${\mu }_{\mathbf{h}}^d\in {\mathbb{C}}^{M{L}_d\times 1}$ and ${\sum }_{\mathbf{h}}^d\in {\mathbb{C}}^{M{L}_d\times M{L}_d}$ denote the corresponding dth cluster in $\mathbf{h}\in {\mathbb{C}}^{ML\times 1}$, ${\mu }_{\mathbf{h}}\in {\mathbb{C}}^{ML\times 1}$ and ${\sum }_{\mathbf{h}}\in {\mathbb{C}}^{ML\times ML}$, respectively. ${\mu }_{\mathbf{h}}$ and ${\sum }_{\mathbf{h}}$ are evaluated according to (19) and (20), given the estimated hyperparameters ${\mathbf{\theta }}^{\left(old\right)}$ and the noise variance $\lambda$. As shown in (14) and (16), ${\mathsf{\Gamma }}_d$ and ${\mathbf{B}}_d$ are determined by the hyperparameters ${\gamma }_d$ and ${\beta }_d$, respectively.

The partial derivative of (24) with respect to ${\gamma }_{d,i}$$(d=1,2,\dots ,C;i=1,2\dots ,L_d)$ is given by

$$\frac{\partial Q}{\partial {\gamma }_{d,i}}=-\frac{M}{2{\gamma }_{d,i}}+\frac{1}{2{\gamma }_{d,i}^2}Tr\left[{\mathbf{B}}_d^{-1}\left({\sum }_{\mathbf{h}}^{d,i}+{\mu }_{\mathbf{h}}^{d,i}{\left({\mu }_{\mathbf{h}}^{d,i}\right)}^H\right)\right],$$

(25)

where ${\mu }_{\mathbf{h}}^{d,i}\in {\mathbb{C}}^{M\times 1}$ denotes the corresponding ith path in ${\mu }_{\mathbf{h}}^d$, and ${\sum }_{\mathbf{h}}^{d,i}\in {\mathbb{C}}^{M\times M}$ is the corresponding principal diagonal block for the ith path in ${\sum }_{\mathbf{h}}^d$. So the learning rule for ${\gamma }_{d,i}$ can be given by

$${\gamma }_{d,i}=\frac{1}{M}Tr\left[{\mathbf{B}}_d^{-1}\left({\sum }_{\mathbf{h}}^{d,i}+{\mu }_{\mathbf{h}}^{d,i}{\left({\mu }_{\mathbf{h}}^{d,i}\right)}^H\right)\right].$$

(26)

The gradient of (24) over ${\mathbf{B}}_d$$(d=1,2,\dots ,C)$ is:

$$\frac{\partial Q}{\partial {\mathbf{B}}_d}=-\frac{L_d}{2}{\mathbf{B}}_d^{-1}+\frac{1}{2}\sum _{i=1}^{L_d}\frac{1}{{\gamma }_{d,i}}{\mathbf{B}}_d^{-1}\left({\sum }_{\mathbf{h}}^{d,i}+{\mu }_{\mathbf{h}}^{d,i}{\left({\mu }_{\mathbf{h}}^{d,i}\right)}^H\right){\mathbf{B}}_d^{-1}.$$

(27)

Thus, the learning rule for ${\mathbf{B}}_d$ can be derived as

$${\mathbf{B}}_d=\frac{1}{L_d}\sum _{i=1}^{L_d}\frac{{\sum }_{\mathbf{h}}^{d,i}+{\mu }_{\mathbf{h}}^{d,i}{\left({\mu }_{\mathbf{h}}^{d,i}\right)}^H}{{\gamma }_{d,i}}.$$

(28)

here, we constrain ${\mathbf{B}}_d$ to have the form as shown in (16), and thus ${\beta }_d$ can be empirically calculated as ${\beta }_d\triangleq \frac{{\alpha }_1}{{\alpha }_0}$ where ${\alpha }_0$ (resp. ${\alpha }_1$) is the average of the elements along the main diagonal (resp. the main sub-diagonal) of the matrix ${\mathbf{B}}_d$ in (28).

4.3. Partition-Based Clustered-Sparse Bayesian Learning Algorithm

When the cluster partition information is known, the hyperparameters $\mathbf{\theta }$, the posterior mean ${\mu }_{\mathbf{h}}$ and the posterior covariance ${\sum }_{\mathbf{h}}$ can be iteratively obtained through the EM method. In other words, cluster partition and cluster distribution jointly determine the clustered-sparse channel distribution. However, it is usually difficult to obtain the cluster partition in advance. As shown in Figure 1, several significant clusters with high energy in CIRs are scattered over a long time spread, and usually do not overlap each other. As shown in (28), the boundary information of these clusters is important for calculating ${\mathbf{B}}_d$, because different clusters have different coherence times (i.e., different AR coefficients). In particular, some weak clusters that can be ignored will cause the hyperparameter estimation error of significant clusters when these clusters are not distinguished. Therefore, the PB-CSBL algorithm is proposed to update cluster partition by pruning the paths with very small energy (i.e., ${\gamma }_{d,i}\to 0$). As shown in Algorithm 1, the PB-CSBL algorithm mainly includes three sub-parts: Channel Estimation, Cluster Partition, Cluster Evolution, and they will loop until the iteration stop condition is reached. In the sub-part of Cluster Partition, only the dense paths with high energy are retained to form clusters and the paths with negligible energy are pruned. For a cluster with some pruned paths, it can further split into more clusters according to the positions of the pruned paths. To reduce the impact of an inaccurate clustering structure at the beginning of the iteration, the maximum possible number of discrete paths in one cluster is given. In the sub-part of Cluster Evolution, the derived cluster distribution learning rules are included. A good initial estimation of the unknown hyperparameters is significant for achieving the global maximum instead of a local maximum. For practical implementation, it was found that an initial estimation given by ${\mathsf{\Gamma }}_d={\mathbf{I}}_{L_d}$ and ${\mathbf{B}}_d={\mathbf{I}}_M$ was sufficient for the proposed PB-CSBL algorithm.

Algorithm 1: Partition-Based Clustered-Sparse Bayesian Learning Algorithm

Input: the received pilot signal ${\mathbf{y}}_{\mathbf{p}}$; the dictionary matrix ${\mathbf{U}}_{\mathbf{P}}$; the noise variance $\lambda$;

the length of discrete paths L; the maximum number of iterations $r_{max}$;

the maximum number of discrete paths in one cluster $L{C}_{max}$; the threshold for

prunning small hyperparameters ${\gamma }_{th}$; the threshold to stop the whole

algorithm $ϵ$.

Initialize: the list of path power ${\gamma }_{list}={\mathbf{1}}_L$; the number of clusters $C=1$;

the list of cluster structure $\mathbf{CLS}$ : ${\mathsf{\Gamma }}_1={\mathbf{I}}_L$, ${\mathbf{B}}_1={\mathbf{I}}_M$ and the delay range of

the 1st cluster $R_d=[0,L-1]$; the iteration counter $r=0$.

Channel Estimation:

1: ${\sum }_{\mathbf{0}}\leftarrow diag\left({\mathsf{\Gamma }}_1\otimes {\mathbf{B}}_1,{\mathsf{\Gamma }}_2\otimes {\mathbf{B}}_2,\dots ,{\mathsf{\Gamma }}_C\otimes {\mathbf{B}}_C\right)$.

2: ${\sum }_{\mathbf{h}}={\sum }_0-{\sum }_0{\mathbf{U}}_{\mathbf{P}}^H{\left(\lambda \mathbf{I}+{\mathbf{U}}_{\mathbf{P}}{\sum }_0{\mathbf{U}}_{\mathbf{P}}^H\right)}^{-1}{\mathbf{U}}_{\mathbf{P}}{\sum }_0$.

3: ${\mu }_{\mathbf{h}}\leftarrow {\lambda }^{-1}{\sum }_{\mathbf{h}}{\mathbf{U}}_{\mathbf{P}}^H{\mathbf{y}}_{\mathbf{p}}$.

Cluster Evolution:

4: for $d=1,2,\dots ,C$do

5:  for  $i=1,2,\dots ,L_d$do

6:   ${\gamma }_{d,i}\leftarrow \frac{1}{M}Tr\left[{\mathbf{B}}_d^{-1}\left({\sum }_{\mathbf{h}}^{d,i}+{\mu }_{\mathbf{h}}^{d,i}{\left({\mu }_{\mathbf{h}}^{d,i}\right)}^H\right)\right]$, and update ${\gamma }_{list}$ with ${\gamma }_{d,i}$.

7:  end for

8:  $L_d^{\prime }\leftarrow min(L_d,L{C}_{max})$.

9:  ${\mathbf{B}}_d\leftarrow \frac{1}{L_d^{\prime }}{\displaystyle \sum _{i=1}^{L_d^{\prime }}}\frac{{\sum }_{\mathbf{h}}^{d,i}+{\mu }_{\mathbf{h}}^{d,i}{\left({\mu }_{\mathbf{h}}^{d,i}\right)}^H}{{\gamma }_{d,i}}$ using the $L_d^{\prime }$ most significant continuous paths.

10:  ${\beta }_d\leftarrow \frac{{\alpha }_1}{{\alpha }_0}$, where ${\alpha }_1$ and ${\alpha }_0$ can be obtained through ${\mathbf{B}}_d$.

11: end for

Cluster Partition:

12: ${\mathbf{p}}^{\left(0\right)}\leftarrow$ FindIndex$(0<{\gamma }_{list}<{\gamma }_{th})$.

13: if ${\mathbf{p}}^{\left(0\right)}$ is not empty then

14:  ${\gamma }_{list}^{\left(old\right)}\leftarrow {\gamma }_{list}$ and ${\gamma }_{list}\left[{\mathbf{p}}^{\left(0\right)}\right]\leftarrow \mathbf{0}$.

15:  $\mathbf{CLS}\leftarrow$ Split($\mathbf{CLS}$, ${\mathbf{p}}^{\left(0\right)}$), where $\mathbf{CLS}$ is splitted into C clusters

according to ${\mathbf{p}}^{\left(0\right)}$.

16:  ${\mathbf{p}}^{\left(1\right)}\leftarrow$ FindIndex$({\gamma }_{list}\ge {\gamma }_{th})$ and ${\mathbf{U}}_{\mathbf{P}}\leftarrow {\mathbf{U}}_{\mathbf{P}}[:,{\mathbf{p}}^{\left(1\right)}\otimes (1:M)]$.

17: end if

18: for $d=1,2,\dots ,C$do

19:  Update $R_d$, ${\mathsf{\Gamma }}_d$ and ${\mathbf{B}}_d$ in $\mathbf{CLS}$ according to ${\mathbf{p}}^{\left(0\right)}$, ${\gamma }_{list}$ and ${\beta }_d$, respectively.

20: end for

Check stopping conditions:

21: $r\leftarrow r+1$.

22: return the sub-part of $\mathbf{Channel}\mathbf{Estimation}$ until $r\ge r_{max}$ or

$\parallel {\gamma }_{list}-{\gamma }_{list}^{\left(old\right)}{\parallel }_2^2<ϵ$.

Output: the pruned channel ${\mu }_{\mathbf{h}}$ with the covariance matrix ${\sum }_{\mathbf{h}}$ and the

cluster list $\mathbf{CLS}$.

4.4. Complexity and Performance Analysis

Considering the excessive delay spread of the UWA channel and the purpose of saving pilot overhead, we assumed $K_d\ge L\ge K_p\ge M$. The time complexity of the proposed PB-CSBL algorithm is compared with the least squares (LS) [39], OMP [25], SOMP [27] and TMSBL [32]. In order to avoid involving matrix inversion, pilot subcarriers are equispaced in the LS method. For fair comparison, the residual Doppler expansion is ignored when evaluating the time complexity of the OMP method. The comparison of time complexity is shown in

Table 1. The comparison of time complexity.

LS	OMP	SOMP	TMSBL	PB-CSBL
$\mathcal{O}\left(L{K}_{p}M\right)$	$\mathcal{O}\left(L^2{K}_{p}M\right)$	$\mathcal{O}\left(L^2{K}_{p}M\right)$	$\mathcal{O}\left(L{K}_p^2\overline{N}\right)$	$\mathcal{O}\left(L{K}_p^2{M}^3\overline{N}\right)$

, where $\overline{N}$ denotes the average iteration number. It can be seen that the time complexity of the PB-CSBL algorithm is higher than that of the LS, OMP and SOMP. This is due to the fact that the temporal correlation of CIRs is ignored in the LS and OMP. Although the stability of path delays is exploited in the SOMP, the correlation of path gains is discarded. Therefore, LS, OMP and SOMP may exhibit worse channel estimation performance compared to the PB-CSBL algorithm. As for the TMSBL with $\mathcal{O}\left(L{K}_p^2\right)$ per iteration, it has lower time complexity than the proposed PB-CSBL with $\mathcal{O}\left(L{K}_p^2{M}^3\right)$ per iteration. This is mainly because it adopts the measure in [41], assuming that there is no temporal correlation among consecutive CIRs or no additive noise in received signal, but this may bring some performance loss, especially for slow-varying or high-noise channels. Therefore, when M is small, we do not adopt the simplification in [41]. The computational complexity of the PB-CSBL algorithm can be further reduced by implementing a first-order algorithm in E-step [34], but the implementation details are omitted here.

According to (20), the mean square error (MSE) bound of the sparse vector $\mathbf{h}$ can be expressed as

$$E\left[\parallel \mathbf{h}-\widehat{\mathbf{h}}{\parallel }_2^2\right]\ge Tr\left({\sum }_{\mathbf{h}}\right)=Tr\left[{\sum }_0-{\sum }_0{\mathbf{U}}_{\mathbf{P}}^H{\left(\lambda \mathbf{I}+{\mathbf{U}}_{\mathbf{P}}{\sum }_0{\mathbf{U}}_{\mathbf{P}}^H\right)}^{-1}{\mathbf{U}}_{\mathbf{P}}{\sum }_0\right],$$

(29)

where $\widehat{\mathbf{h}}$ is the estimate of $\mathbf{h}$, and its estimation accuracy is related to the hyperparameter vector $\mathbf{\theta }$ obtained by maximizing the marginal likelihood function. Thus, the lower bound of MSE is reached (i.e., equality holds) only when the perfect hyperparameters are obtained.

Assume that all L paths are in one cluster (i.e., $C=1$) and consider two special cases of $\mathbf{\theta }$:

(1)

The temporal correlation coefficient $\beta =0$.

In this case, there is no temporal correlation among CIRs. It is assumed that each transmitted symbol is normalized to unit power. Then, the MSE bound of ${\mathbf{h}}^{\left(m\right)}$ can be expressed by

$$E\left[\parallel {\mathbf{h}}^{\left(m\right)}-{\widehat{\mathbf{h}}}^{\left(m\right)}{\parallel }_2^2\right]\ge Tr\left[{({\mathsf{\Gamma }}^{-1}+{\lambda }^{-1}{\mathsf{\Phi }}_{\mathbf{P}}^H{\mathsf{\Phi }}_{\mathbf{P}})}^{-1}\right]=\sum _{l=1}^L{\left(\frac{K_p}{\lambda }+\frac{1}{{\gamma }_l}\right)}^{-1},$$

(30)

where $\mathsf{\Gamma }=diag\left(\gamma \right)$, $\gamma =[{\gamma }_1,{\gamma }_2,\dots ,{\gamma }_L]$ and ${\gamma }_l$ controls the variance of the lth channel coefficient.

(2)

The temporal correlation coefficient $\beta =1$

In this case, the CIRs are time-invariant across M OFDM blocks. Thus, the additive noise variance $\lambda$ can be reduced to $\lambda /M$, and the MSE bound of ${\mathbf{h}}^{\left(m\right)}$ can be expressed by

$$E\left[\parallel {\mathbf{h}}^{\left(m\right)}-{\widehat{\mathbf{h}}}^{\left(m\right)}{\parallel }_2^2\right]\ge \sum _{l=1}^L{\left(\frac{M{K}_p}{\lambda }+\frac{1}{{\gamma }_l}\right)}^{-1}.$$

(31)

As shown in the case of the high temporal correlation (i.e., $\beta =1$), the anti-noise performance of channel estimation can be improved. As for the low temporal correlation among CIRs (i.e., $\beta =0$), the anti-noise performance is not significantly improved. However, with more reference information, the better hyperparameter vector $\gamma$ can be obtained. Therefore, when multiple OFDM blocks are utilized simultaneously, higher channel estimation accuracy can be achieved.

5. Noise-Corrected Clustered-Sparse Channel Estimation

In this section, the parametric form of noise-resistant Bayesian channel estimation is given first. Then, the data detection and noise measurement are derived. Then, the noise-corrected clustered-sparse channel estimation algorithm is presented. At last, the complexity and performance analysis of this algorithm is given.

5.1. Noise-Resistant Bayesian Channel Estimation

As the number of pilots $K_p$ decreases, the precision of channel estimation using the PB-CSBL algorithm will rapidly deteriorate. This is primarily because the algorithm does not utilize all the information available. Specifically, although both pilot and data symbols are transmitted, only the $K_p\times M$ matrix ${\mathbf{Y}}_{\mathbf{P}}$ corresponding to the pilot subcarriers in (8) is used for estimating the CIRs. For the matrix $\mathbf{Y}$ in (7), the remaining $K_d\times M$ observations are typically not used as they contain the unknown data symbols. If the unknown data symbols are perfectly obtained, using the Bayes rule, the posterior mean of CIRs based on the prior probability in (17) and the channel measurement equation in (10) can be written by

$${\overline{\mu }}_{\mathbf{h}}={\sum }_0{\mathbf{U}}^H{\left(\mathsf{\Lambda }+\mathbf{U}{\sum }_0{\mathbf{U}}^H\right)}^{-1}\mathbf{y},$$

(32)

with the posterior covariance matrix:

$${{\displaystyle \overline{\sum }}}_{\mathbf{h}}={\sum }_0-{\sum }_0{\mathbf{U}}^H{\left(\mathsf{\Lambda }+\mathbf{U}{\sum }_0{\mathbf{U}}^H\right)}^{-1}\mathbf{U}{\sum }_0.$$

(33)

Moreover, AWGN is assumed to simplify implementation complexity in the PB-CSBL algorithm, as the noise power distribution $\mathbf{\lambda }=[{\lambda }_{-K/2},{\lambda }_{-K/2+1},\dots ,{\lambda }_{K/2-1}]$ depends on the observation of all subcarriers. However, in a complex UWA environment, the noise power may vary significantly over the signal bandwidth. Using decided data symbols and corrected colored noise, the high estimation accuracy can be achieved by the MAP estimate in (32).

5.2. Data Detection and Noise Measurement

To obtain the MAP estimate in (32), the hyperparameters related to the cluster distribution, $\mathbf{\theta }\triangleq \left[{\gamma }_1,{\gamma }_2,\dots ,{\gamma }_C,{\beta }_1,{\beta }_2,\dots ,{\beta }_C\right]$, can be iteratively updated, as shown in the proposed PB-CSBL algorithm, however, the detection of data signal ${\mathbf{x}}_{\mathbf{d}}^{\left(m\right)}$ and the noise distribution $\mathbf{\lambda }$ are still unresolved. Using the notations ${\mathbf{x}}_{\mathbf{d}}\triangleq \left[{\left({\mathbf{x}}_{\mathbf{d}}^{\left(1\right)}\right)}^T,{\left({\mathbf{x}}_{\mathbf{d}}^{\left(2\right)}\right)}^T,\dots ,{\left({\mathbf{x}}_{\mathbf{d}}^{\left(M\right)}\right)}^T\right]$, $\mathbf{\phi }\triangleq \left[\lambda ,{\mathbf{x}}_{\mathbf{d}}\right]$, and $\mathbf{\zeta }\triangleq \left[\mathbf{\theta },\mathbf{\phi }\right]$, the unknown hyperparameters set $\mathbf{\zeta }$ can be estimated by maximizing the marginal likelihood function $p(\mathbf{y};\mathbf{\zeta })$. This is equivalent to minimizing $-logp(\mathbf{y};\mathbf{\zeta })$, giving the effective cost function:

$$\mathcal{L}\left(\mathbf{\zeta }\right)={\mathbf{y}}^H{\sum }_{\mathbf{y}}^{-1}\left(\mathbf{\zeta }\right)\mathbf{y}+log\left|{\sum }_{\mathbf{y}}\left(\mathbf{\zeta }\right)\right|,$$

(34)

where ${\sum }_{\mathbf{y}}\left(\mathbf{\zeta }\right)=\mathsf{\Lambda }\left(\lambda \right)+\mathbf{U}\left({\mathbf{x}}_{\mathbf{d}}\right){\sum }_0\left(\mathbf{\theta }\right){\mathbf{U}}^H\left({\mathbf{x}}_{\mathbf{d}}\right)$. Similar to the cluster distribution learning, the expectation maximization (EM) method [41] can be employed to solve (34) iteratively. While $\mathbf{h}$ is considered as the nuisance variable, $\mathbf{\zeta }$ can be estimated by maximizing:

$$Q\left(\mathbf{\zeta }\right|{\mathbf{\zeta }}^{\left(old\right)})=E_{\mathbf{h}|\mathbf{y};{\mathbf{\zeta }}^{\left(old\right)}}\left[logp(\mathbf{y},\mathbf{h};\mathbf{\zeta })\right]=E_{\mathbf{h}|\mathbf{y};{\mathbf{\zeta }}^{\left(old\right)}}\left[logp\left(\mathbf{y}\right|\mathbf{h};\mathbf{\phi })\right]+E_{\mathbf{h}|\mathbf{y};{\mathbf{\zeta }}^{\left(old\right)}}\left[logp(\mathbf{h};\mathbf{\theta })\right],$$

(35)

where ${\mathbf{\zeta }}^{\left(old\right)}$ denotes the estimated hyperparameters in the previous iteration. As shown in (35), the joint maximization simplifies into two independent maximizations over $\mathbf{\theta }$ and $\mathbf{\phi }$. On optimizing the latter function in (35) with respect to $\mathbf{\theta }$, one obtains the learning rules about $\mathbf{\theta }$ as (26) and (28). Maximizing the former over ${\mathbf{x}}_{\mathbf{d}}^{\left(m\right)}$ yields the minimum mean squared error (MMSE) equalizer:

$${\widehat{\mathbf{x}}}_{\mathbf{d}}^{\left(m\right)}={\left[diag\left({\lambda }_{\mathbf{d}}\right){\mathbf{C}}_{\mathbf{X}}^{-1}+Prdiag\left({\mathbf{F}}_{\mathbf{D}}{\Xi }_{\mathbf{h}}^{\left(m\right)}{\mathbf{F}}_{\mathbf{D}}^H\right)\right]}^{-1}diag\left[{\left({\overline{\mu }}_{\mathbf{h}}^{\left(m\right)}\right)}^H{\mathbf{F}}_{\mathbf{D}}^H\right]{\mathbf{y}}_{\mathbf{d}}^{\left(m\right)},$$

(36)

where ${\mathbf{y}}_{\mathbf{d}}^{\left(m\right)}$ and ${\lambda }_{\mathbf{d}}$ are the subvectors of ${\mathbf{y}}^{\left(m\right)}$ and $\lambda$ at the location of data subcarriers, respectively. ${\mathbf{C}}_{\mathbf{X}}$ is the covariance matrix of ${\mathbf{x}}_{\mathbf{d}}^{\left(m\right)}$, ${\mathbf{F}}_{\mathbf{D}}$ is the $K_d\times L$ DFT matrix at the location of data subcarriers. ${\Xi }_{\mathbf{h}}^{\left(m\right)}$ is defined as ${\Xi }_{\mathbf{h}}^{\left(m\right)}\triangleq {\overline{\mu }}_{\mathbf{h}}^{\left(m\right)}{\left({\overline{\mu }}_{\mathbf{h}}^{\left(m\right)}\right)}^H+{{\displaystyle \overline{\sum }}}_{\mathbf{h}}^{\left(m\right)}$, where ${\overline{\mu }}_{\mathbf{h}}^{\left(m\right)}\in {\mathbb{C}}^{L\times 1}$ and ${{\displaystyle \overline{\sum }}}_{\mathbf{h}}^{\left(m\right)}\in {\mathbb{C}}^{L\times L}$ denote the posterior mean and covariance matrix of the mth OFDM block, respectively. The decisions on the data symbols, ${\mathbf{x}}_{\mathbf{d}}^{\left(m\right)}$, are obtained by mapping the soft decisions, ${\widehat{\mathbf{x}}}_{\mathbf{d}}^{\left(m\right)}$, to the closest point on the constellation, meaning that:

$${\mathbf{x}}_{\mathbf{d}}^{\left(m\right)}=Quantization\left({\widehat{\mathbf{x}}}_{\mathbf{d}}^{\left(m\right)}\right).$$

(37)

Although the noise distribution is non-uniform, the noise powers of adjacent subcarriers are usually close (especially when the number of subcarriers is large). Then, we divide all subcarriers into non-overlapping Q groups, so that $K/Q$ subcarriers can be grouped into the same one. Thus, the average noise power of the qth group can be calculated as

$${\lambda }_g^{\left(q\right)}=\frac{Q}{KM}Tr\left[{\mathbf{y}}^{\left(q\right)}{\left({\mathbf{y}}^{\left(q\right)}\right)}^H-{\mathbf{U}}^{\left(q\right)}{\Xi }_{\mathbf{h}}{\left({\mathbf{U}}^{\left(q\right)}\right)}^H\right],$$

(38)

where ${\Xi }_{\mathbf{h}}$ is defined as ${\Xi }_{\mathbf{h}}\triangleq {\overline{\mu }}_{\mathbf{h}}{\left({\overline{\mu }}_{\mathbf{h}}\right)}^H+{{\displaystyle \overline{\sum }}}_{\mathbf{h}}$. $\mathbf{U}\in {\mathbb{C}}^{KM\times LM}$ can be calculated according to (11), as unknown data symbols are estimated by (36) and (37). ${\mathbf{U}}^{\left(q\right)}\in {\mathbb{C}}^{\frac{KM}{Q}\times LM}$ is the submatrix of $\mathbf{U}$ and denotes the qth group. ${\mathbf{y}}^{\left(q\right)}$ denotes the qth group of $\mathbf{y}$. After ${\lambda }_g^{\left(1\right)},{\lambda }_g^{\left(2\right)},\dots ,{\lambda }_g^{\left(Q\right)}$ are calculated according to (38), the non-uniform noise power distribution across all subcarriers, $\mathbf{\lambda }=[{\lambda }_{-K/2},{\lambda }_{-K/2+1},\dots ,{\lambda }_{K/2-1}]$, can be obtained by upsampling the Q average noise powers.

5.3. Noise-Corrected Clustered-Sparse Channel Estimation Algorithm

Based on the noise-resistant Bayesian channel estimation, noise measurement and data detection, we can develop the PB-CSBL algorithm into the NC-CSCE Algorithm. As shown in Algorithm 2, the sub-parts of Cluster Partition and Cluster Evolution have the similar process to the PB-CSBL algorithm. With the aid of Data Detection and Noise Measurement, the estimation accuracy in Channel Estimation is improved. Conversely, the improvement of channel estimation also ensures the high accuracy of data detection and noise measurement results. To achieve the global maximum instead of a local maximum, the NC-CSCE algorithm requires a good initial estimate with respect to the unknown cluster distribution, data signal and noise distribution. However, it is obviously difficult to initialize so many hyperparameters at the same time. Therefore, the PB-CSBL algorithm is utilized to initialize ${\overline{\mu }}_{\mathbf{h}}$ and ${{\displaystyle \overline{\sum }}}_{\mathbf{h}}$ instead of initializing these hyperparameters. Once ${\overline{\mu }}_{\mathbf{h}}$ and ${{\displaystyle \overline{\sum }}}_{\mathbf{h}}$ are initialized, these hyperparameters can be calculated as shown in Algorithm 2.

5.4. Complexity and Performance Analysis

The NC-CSCE algorithm utilizes all available subcarriers to achieve better channel estimation performance. Therefore, compared with the proposed PB-CSBL algorithm with $\mathcal{O}\left(L{K}_p^2{M}^3\right)$ per iteration, the NC-CSCE algorithm has higher time complexity, i.e., $\mathcal{O}\left({(K_d+K_p)}^3{M}^3\right)$ per iteration. Similar to the PB-CSBL algorithm, the complexity of the NC-CSCE algorithm can be further simplified by using the techniques in [34,41], however, these simplifications may lead to the performance degradation of channel estimation because of the errors introduced by the approximations. The details of the simplifications are omitted here.

Algorithm 2: Noise-Corrected Clustered-Sparse Channel Estimation Algorithm

Input: the received signal $\mathbf{y}$; the noise variance $\lambda$; the length of discreted paths L;

the maximum number of iterations $r_{max}$; the maximun number of discreted

paths in one cluster $L{C}_{max}$; the threshold for prunning small

hyperparameters ${\gamma }_{th}$; the threshold to stop the whole algorithm $ϵ$.

Initialize: The PB-CSBL algorithm is utilized to obtain the posterior mean ${\overline{\mu }}_{\mathbf{h}}\leftarrow {\mu }_{\mathbf{h}}$

and covariance matrix ${{\displaystyle \overline{\sum }}}_{\mathbf{h}}\leftarrow {\sum }_{\mathbf{h}}$; the noise vector $\lambda =\lambda {\mathbf{I}}_K$; the list of path

power ${\gamma }_{list}={\mathbf{1}}_L$; the iteration counter $r=0$.

Cluster Evolution:

1: for $d=1,2,\dots ,C$do

2:  for $i=1,2,\dots ,L_d$do

3:   ${\gamma }_{d,i}\leftarrow \frac{1}{M}Tr\left[{\mathbf{B}}_d^{-1}\left({{\displaystyle \overline{\sum }}}_{\mathbf{h}}^{d,i}+{\overline{\mu }}_{\mathbf{h}}^{d,i}{\left({\overline{\mu }}_{\mathbf{h}}^{d,i}\right)}^H\right)\right]$, and update ${\gamma }_{list}$ with ${\gamma }_{d,i}$.

4:  end for

5:  $L_d^{\prime }\leftarrow min(L_d,L{C}_{max})$.

6:  ${\mathbf{B}}_d\leftarrow \frac{1}{L_d^{\prime }}\sum _{i=1}^{L_d^{\prime }}\frac{{{\displaystyle \overline{\sum }}}_{\mathbf{h}}^{d,i}+{\overline{\mu }}_{\mathbf{h}}^{d,i}{\left({\overline{\mu }}_{\mathbf{h}}^{d,i}\right)}^H}{{\gamma }_{d,i}}$ using the $L_d^{\prime }$ most significant continuous paths.

7:  ${\beta }_d\leftarrow \frac{{\alpha }_1}{{\alpha }_0}$, where ${\alpha }_1$ and ${\alpha }_0$ can be obtained through ${\mathbf{B}}_d$.

8: end for

Cluster Partition:

9: ${\mathbf{p}}^{\left(0\right)}\leftarrow$ FindIndex$(0<{\gamma }_{list}<{\gamma }_{th})$.

10: if${\mathbf{p}}^{\left(0\right)}$ is not empty then

11:  ${\gamma }_{list}^{\left(old\right)}\leftarrow {\gamma }_{list}$ and ${\gamma }_{list}\left[{\mathbf{p}}^{\left(0\right)}\right]\leftarrow \mathbf{0}$.

12:  $\mathbf{CLS}\leftarrow$ Split($\mathbf{CLS}$, ${\mathbf{p}}^{\left(0\right)}$), where $\mathbf{CLS}$ is splitted into C clusters according

to ${\mathbf{p}}^{\left(0\right)}$.

13: end if

14: for $d=1,2,\dots ,C$do

15:  Update $R_d$, ${\mathsf{\Gamma }}_d$ and ${\mathbf{B}}_d$ in $\mathbf{CLS}$ according to ${\mathbf{p}}^{\left(0\right)}$, ${\gamma }_{list}$ and ${\beta }_d$, respectively.

16: end for

Data Detection:

17: Refer to (36) and (37).

Noise Measurement:

18: Update the noise vector $\lambda$ according to (38).

Channel Estimation:

19: ${\sum }_{\mathbf{0}}\leftarrow diag\left({\mathsf{\Gamma }}_1\otimes {\mathbf{B}}_1,{\mathsf{\Gamma }}_2\otimes {\mathbf{B}}_2,\dots ,{\mathsf{\Gamma }}_C\otimes {\mathbf{B}}_C\right)$.

20: Obtain $\mathbf{U}$ according to (11).

21: ${\mathbf{p}}^{\left(1\right)}\leftarrow$ FindIndex$({\gamma }_{list}\ge {\gamma }_{th})$ and $\mathbf{U}\leftarrow \mathbf{U}[:,{\mathbf{p}}^{\left(1\right)}\otimes (1:M)]$.

22: Refer to (32) and (33).

Check stopping conditions:

23: $r\leftarrow r+1$.

24: return the sub-part of $\mathbf{Cluster}\mathbf{Evolution}$ until $r\ge r_{max}$ or $\parallel {\gamma }_{list}-{\gamma }_{list}^{\left(old\right)}{\parallel }_2^2<ϵ$.

Output: the pruned channel ${\overline{\mu }}_{\mathbf{h}}$ with the covariance matrix ${{\displaystyle \overline{\sum }}}_{\mathbf{h}}$, the cluster list $\mathbf{CLS}$,

the noise vector $\lambda$, and the transmitted data symbols ${\mathbf{x}}_{\mathbf{d}}^{\left(m\right)}$ for $m=1,2,\dots ,M$.

Based on the noise-resistant Bayesian channel estimation in the NC-CSCE algorithm, the MSE bound of the sparse vector $\mathbf{h}$ can be expressed as

$$E\left[\parallel \mathbf{h}-\widehat{\overline{\mathbf{h}}}{\parallel }_2^2\right]\ge Tr\left({{\displaystyle \overline{\sum }}}_{\mathbf{h}}\right)=Tr\left[{\sum }_0-{\sum }_0{\mathbf{U}}^H{\left(\mathsf{\Lambda }+\mathbf{U}{\sum }_0{\mathbf{U}}^H\right)}^{-1}\mathbf{U}{\sum }_0\right],$$

(39)

where $\widehat{\overline{\mathbf{h}}}$ is the estimate of $\mathbf{h}$. The MSE bound is related to the hyperparameters of the cluster distribution, data symbols and noise distribution. Only when the perfect hyperparameters are obtained is the lower bound of MSE reached. Similarly to analyzing the PB-CSBL algorithm, we assume that all L paths are in one cluster (i.e., $C=1$) and consider the two special cases, namely the temporal correlation coefficient $\beta =0$ and $\beta =1$. For ease of analysis, the noise power keeps constant $\overline{\lambda }$ over all subcarriers. Furthermore, each transmitted symbol is normalized to unit power. Then, the MSE bounds of ${\mathbf{h}}^{\left(m\right)}$ can be evaluated by

$$E\left[\parallel {\mathbf{h}}^{\left(m\right)}-{\widehat{\overline{\mathbf{h}}}}^{\left(m\right)}{\parallel }_2^2\right]\ge \sum _{l=1}^L{\left(\frac{K_p+K_d}{\overline{\lambda }}+\frac{1}{{\gamma }_l}\right)}^{-1},$$

(40)

$$E\left[\parallel {\mathbf{h}}^{\left(m\right)}-{\widehat{\overline{\mathbf{h}}}}^{\left(m\right)}{\parallel }_2^2\right]\ge \sum _{l=1}^L{\left(\frac{M(K_p+K_d)}{\overline{\lambda }}+\frac{1}{{\gamma }_l}\right)}^{-1},$$

(41)

respectively. In addition to what has been discussed when analyzing the PB-CSBL algorithm, Equations (40) and (41) also prove that the NC-CSCE algorithm utilizing pilot symbols and unknown data symbols can achieve a smaller MSE compared with the PB-CSBL algorithm.

6. Evaluation and Result Analysis

In this section, to evaluate the performance of the proposed algorithms, numerical simulations and a lake trial were conducted. The parameter settings of CP-OFDM in the simulations and lake trial are shown in

Table 2. UWA CP-OFDM settings.

Parameters	Notations	Values of the Simulations	Values of the Lake Trial
Bandwidth	B	1.6 kHz	5 kHz
Carrier frequency	$f_c$	2.5 kHz	20 kHz
Sampling frequency	$f_s$	12.5 kHz	100 kHz
Number of subcarriers	K	256	256
Number of data subcarriers	$K_d$	218	193
Number of pilot subcarriers	$K_p$	14	32
Number of null subcarriers	$K_n$	24	31
Symbol duration without CP	T	160 ms	51.2 ms
CP length	$T_{cp}$	10 ms	25.6 ms
Blocks in one frame	M	4	5

. The bandwidth, carrier frequency and sampling frequency of the lake trial are closely related to the transducers used in the trial, and these values are set smaller in the simulation system to increase the simulation speed. Because of the high cost of the lake trial, the CP length, the number of pilot subcarriers and the number of null subcarriers are designed to be larger than the simulation system to ensure the availability of trial data. In the simulation system (the specifications of our computer are as follows: Intel(R) Core(TM) i5-6400 CPU 2.7 GHz (4 cores), 16 GB RAM, 1 TB memory with Windows and MATLAB installed), it is assumed that the UWA channel has 12 randomly generated paths, and that these paths remain fixed within one OFDM frame. The maximum path delay is less than CP length, and the amplitudes of paths are Rayleigh distributed with the average power decreasing exponentially with delay. The Doppler scale factor $a\in (0,{10}^{-3})$ is randomly chosen. Moreover, there is no channel code module in the simulation system.

6.1. Simulation Results

Firstly, the proposed PB-CSBL and NC-CSCE algorithms were compared with the LS, OMP, SOMP and TMSBL for channel estimation. All paths are distributed in three non-overlapping clusters, and these clusters are set with three different temporal correlation coefficients ranges over multiple OFDM blocks, which are $[0.1,0.3]$, $[0.5,0.7]$, $[0.8,1.0]$. The additive noise is assumed to be white Gaussian noise, and thus all subcarriers are divided into the same group (i.e., $Q=1$) for the proposed NC-CSCE algorithm. The transmitted data bits are mapped to QPSK or 16-QAM constellation symbols. The performance metrics used in the numerical simulations are the MSE of CIRs and BER of data bits, which are verified in 500 simulations for each value of SNR.

Figure 3a,b show the MSE and BER performance of different channel estimation methods, respectively. LS performs worst, because the channel sparsity is not considered for this underdetermined system. Based on the joint estimation of multiple OFDM blocks, the SOMP makes full use of the strong temporal correlation of path delays and achieves better performance than the OMP. The delay and correlation of path gains are analyzed in the TMSBL, which makes it better than the SOMP. However, the proposed PB-CSBL performs even better, as it can learn different temporal correlation coefficients and reduce the interference among clusters. The proposed NC-CSCE algorithm utilizes all received pilot and data symbols for the channel estimation, which achieves the lowest BER close to the perfect CSI. When 16-QAM mapping is adopted, Figure 4 proves that the proposed algorithms still maintain good estimation performance. However, in the low SNR region from 0 to 10 dB, the NC-CSCE algorithm does not perform well. This is mainly because high-order modulation at low SNR is prone to symbol decision errors, which may be exacerbated in the iterative process of the NC-CSCE algorithm.

Then, to verify the estimation performance of the proposed algorithms under extreme temporal correlation coefficients (i.e., $\beta =0$ and $\beta =1$), the proposed algorithms are compared with their theoretical lower bounds. Multipaths are assumed to be sparse and share the same correlation coefficient. Then, under QPSK mapping and white Gaussian noise, the mean square error performance is presented in Figure 5. When $\beta =0$, the MSE of the proposed algorithms is close to the theoretical lower bound in a high SNR region. A significant gap exists in the low SNR region, especially for the NC-CSCE algorithm. This is caused by imperfect hyperparameter estimation, as the lower bound is only reached when the perfect hyperparameters are obtained. When $\beta =1$, the gap is more obvious in the whole SNR region because of the hyperparameter estimation error. However, the MSE performance with $\beta =1$ is better than the MSE performance with $\beta =0$ in a low SNR region. In short, a larger temporal correlation coefficient means stronger resistance to additive noise, however, the hyperparameter estimation error, such as the error of $\gamma$, is often larger than that of small temporal correlation coefficients in the high SNR region. Figure 6a,b show the MSE of the $\gamma$ estimated by the proposed PB-CSBL and NC-CSCE algorithms, respectively. The MSE is evaluated under different SNRs and temporal correlation coefficients. As can be seen, a larger temporal correlation coefficient achieves a smaller MSE in the low SNR region, while a smaller temporal correlation coefficient reduces the MSE value in the high SNR region.

Finally, to evaluate the estimation performance of the proposed NC-CSCE algorithm under colored noise, we compared the NC-CSCE algorithm with the TMSBL and PB-CSBL algorithms. The group number of subcarriers in NC-CSCE is adjusted to evaluate the channel estimation performance. In the simulation, all paths are distributed in three non-overlapping clusters. The clusters are set with three different temporal correlation coefficients ranges over multiple OFDM blocks, which are $[0.1,0.3]$, $[0.5,0.7]$, $[0.8,1.0]$. All subcarriers are divided into eight groups and each group is set with a different SNR. As shown in Figure 7b, the reference SNR ranges are from 5 to 20 dB. With QPSK mapping, the MSE and BER performance curves are shown in Figure 7a, and the results prove that the proposed NC-CSCE algorithm is superior to the TMSBL and PB-CSBL algorithms. The NC-CSCE-8 denotes the NC-CSCE algorithm with eight groups, and it performs better than the NC-CSCE algorithm with other groups. Therefore, the setting of group number is critical for the NC-CSCE algorithm. In one instance, the actual distribution of colored noise can hardly be reflected with a limited group number; yet in another, a small symbol number will result in low measurement accuracy in one group. Figure 7b also shows the statistical results of the estimated SNRs based on 500 simulations. Similarly, NC-CSCE-8 matches the reference best among all compared methods.

6.2. Lake Experimental Results

To further verify the proposed algorithms, a real lake experiment was performed at Huating Lake, Anhui Province, China, in May 2018. Huating Lake covers an area of about 70 square kilometers, with an average water depth of about 17 m. We hung the transmitter and receiver on the side of a passenger ship and a speedboat, respectively. The passenger ship has a semi-open roof and is about 10 m long, while the speedboat is about half the size of the ship. During the experiment, the transmitter depth was about 10 m and the receiver depth was about 13 m. Affected by wind and current, their relative speed was about 0.5 m/s, and the distance between them was more than 4.5 km. The parameter settings of the experiment are in Table 2. In addition, the linear frequency modulation (LFM) signal was added at the beginning and end of each frame for time synchronization and Doppler estimation. The transmitted data bits were encoded by a $1/2$ rate convolutional code, and the encoded bits were mapped by QPSK constellation. As showed in Figure 1, using the proposed NC-CSCE algorithm, we obtained the estimated CIRs for a 14 frame history with five consecutive OFDM blocks in one frame. The CIRs are estimated after Doppler compensation, and the total delay spread is around 15 ms. From the estimated CIRs, we notice that the paths are distributed in multiple clusters, and the channel shows a clustered-sparse structure.

Figure 8 shows the BER performance after convolutional decoding. The BERs of the proposed two algorithms are still the lowest for most frames, followed by the TMSBL and the SOMP, while the OMP has the worst BER performance. Obviously, the performance gap between proposed algorithms and other methods further verify the superior performance of the proposed algorithms.

7. Conclusions

In this paper, we propose a novel Bayesian learning-based channel estimation architecture to estimate the time-varying multipath channel with colored noise for UWA-OFDM systems. Specifically, a clustered-sparse channel distribution model was constructed to characterize the multipath distribution, and a noise-resistant channel measurement model is constructed to reduce the noise disturbance. To obtain the clustered-sparse distribution, we propose the partition-based clustered-sparse Bayesian learning algorithm. To lessen the effect of colored noise, we proposed a noise-corrected clustered-sparse channel estimation algorithm to improve the estimation performance. Taking advantage of the iterative clustered-sparse distribution learning, symbol decision and noise correction, the accuracy of Bayesian channel estimation can be improved. Experiments proved the effectiveness of the proposed algorithms for channel estimation.

In future work, channel decoding can be incorporated into the Bayesian learning-based channel estimation architecture to explore the possibility of reducing symbol decision errors. In the noise-resistant channel measurement model, the channel gains are approximately constant within one OFDM block, which may be too optimistic in some rapidly time-varying UWA channels. Therefore, some channel estimation models that deal with doubly selective fading channels, such as the popular basis expansion model, can be combined with the clustered-sparse Bayesian learning to improve the estimation performance. At the same time, channel feedback and reconstruction will be analyzed to realize adaptive modulation and coding in UWA-OFDM communication.