A passband transmit cyclic prefix (CP)-OFDM signal in continuous time can be expressed as follows: w ∼ t = 2 Re ∑ k = − K / 2 K / 2 − 1 w k e j 2 π f k t q t , where w k is the transmitted symbol; K is the number of subcarriers in one OFDM block, including data subcarriers K d , pilot subcarriers K p , and null subcarriers K n . f k = f c + k / T is the frequency of the kth subcarriers; and f c is the carrier frequency. Let T denote one OFDM symbol duration, and T c p denote the length of the circular prefix. q t is the pulse-shaping filter: q t = 1 , t ∈ − T cp , T 0 , otherwise .

Figure 1 shows the system model of UWA OFDM communication. Each vertical column in a frame represents an OFDM block; each grid in one vertical column represents each subcarrier frequency; the blue grids indicate the position of the inserted pilots; the white grids indicate the position of the transmitted data. The comb-shaped pilots are integrated into each OFDM block at the same interval in the frequency domain, which can better estimate and compensate for UWA channels in OFDM communication.

In OFDM communication, a frame signal contains multiple OFDM blocks. Each block needs to be estimated and demodulated individually. The cycle prefixes are in the front of the blocks to prevent symbol interference. Pilots are inserted into each block at four equal intervals. The specific frame structure is shown in Figure 2.

We assume that the time-varying underwater acoustic channel is a multipath channel containing M paths; β m and τ m are the gain and delay of the mth path, and a is the Doppler factor. In addition, we assume that these parameters are constant within one CP-OFDM block duration. Then, the channel impulse response can be written as follows: h τ , t = ∑ m = 1 M β m δ τ − τ m − a t .

The received passband signal is as follows: y ∼ t = ∑ m = 1 M β m w ∼ 1 + a t − τ m + z ∼ t , where z ∼ t is the additive noise. The method of CP self-correlators is used to estimate Doppler block by block, and then the received signal is resampled with Doppler factor a . After Doppler estimation and resampling, the Doppler effect of the signal is considered to have been removed. We can model the K × 1 received signal Y for the nth block as follows: Y n = W n F h n + Z n , where W n is the K × K dimensional diagonal matrix consisting of transmitted symbols. Z n is the noise vector. F is the K × M discrete Fourier transform (DFT) matrix. The overall channel is represented as h n = h 1 n , h 2 n , … , h M n , where most of the multipath delay amplitude parameters are 0. One frame consists of N consecutive OFDM blocks, so n ∈ 1 , N . The received model considering only P pilot subcarriers can be written as follows: Y p n = W p n F p h n + Z p n , where Y p n is the K × 1 dimensional pilot reception vector, the diagonal matrix of the transmission pilot is W p n , and Z p n is the noise vector.

The underwater acoustic channel is a typical time-varying sparse channel, with a small number of sparse non-zero paths. The consecutive OFDM symbols have a stable multi-path structure of the channel, which has a temporal correlation for gains. Therefore, we model a sparse channel matrix of N consecutive OFDM symbols as follows: h ¯ ≜ h 1 , … , h n , … , h N , where h n n ∈ 1 , N represents the channel impulse response of the nth OFDM symbol. For each channel vector, the non-zero path shows joint sparsity. The P × M dimensional dictionary matrix is shown as follows: Φ p n = W p n F p .

The consecutive OFDM symbols have the same pilots. Therefore, we can model a joint estimation model for N consecutive OFDM blocks as follows: Y ¯ p n = Φ p h ¯ + Z ¯ p , (1) where Y ¯ p = Y p 1 , … , Y p n , … , Y p N , W ¯ p = W p 1 , … , W p n , … , W p N .

For the estimation problem in Equation 1, we adopt the TMSBL algorithm to jointly estimate the channel matrix h ¯ . First, we model h ¯ i as the ith sparse block of h ¯ , and assume that all sparse signal blocks are independent. The parametric form of the prior of each h ¯ i is as follows: p h ¯ i ; γ i , B i ∼ C N 0 , γ i B i   i = 1 , … , M , where γ i is a non-negative hyperparameter that determines the sparsity of the current block and the channel length is M . B i is a positive definite time correlation matrix, representing the time correlation structure of h ¯ . γ i and B i are unknown hyperparameters that can be estimated in the TMSBL algorithm. We denote Γ as a M × M diagonal matrix with γ . Through the prior distribution of each sparse block, the overall prior distribution of all sparse blocks can be obtained as follows: p h ¯ i ; Γ , B i = ∏ i = 1 M p h ¯ i ; γ i , B i .

According to the prior distribution h ¯ i and likelihood function p Y ¯ p | h ¯ , using Bayesian criterion, we can obtain the posterior probability density of h ¯ as follows: p h ¯ n | Y ¯ p ; Γ = C N μ n , Σ   n ∈ 1 , N .

The covariance and mean are Σ = σ − 2 Φ P H Φ P + Γ r − 1 − 1 , (2) M = μ 1 , … , μ N = σ − 2 Σ Φ P H Y ¯ p . (3)

The μ n is the estimated h n , and the hyperparameter γ i and time correlation matrix B i can be estimated using the EM algorithm, which obtains the parameter by iterative calculation. In the rth iteration, the E step calculates the expectation of variables under conditional probability distribution using Equations 2, 3, and the M-step is expressed via the following updated rule: γ i = 1 M M i B − 1 M i H + Σ i i .

It is worth noting that if each B i is estimated independently, it will lead to overfitting by limited samples and too many parameters; therefore, we use one positive definite matrix B to model all the temporal correlation matrices B i (Zhang and Rao, 2011; Cawley and Nicola, 2007; Guyon et al., 2010): B = ∑ i = 1 M M i H M i γ i + η I N .

The noise variance can be calculated using null subcarriers as follows: σ 2 = E Y ¯ n 2 .

The steps of the TMSBL algorithm are provided in Table 1.

For the estimation problem in Equation 1, we directly list the posterior probability density and likelihood function as follows: p h ¯ n | Y ¯ p ; Γ , B = C N μ n , Σ , B   n ∈ 1 , N , p Y ¯ p | Γ , B = C N 0 , C A , B , where Γ = d i a g − 1 γ i is the M × M diagonal hyperparameter matrix and the covariance and mean are Σ A ‐ 1 = Γ ‐ 1 + σ ‐ 2 Φ P H Φ P and μ = σ − 2 Σ A Φ P H Y , respectively. We define the parameter C A = σ ‐ 2 I M + Φ P H Γ Φ P , where the noise variance is σ 2 . The cost function L of the temporal correlation MMV model can be obtained using the type-2 maximum likelihood estimation method (Tipping, 2001): L = M l o g B + N l o g C A + T r B − 1 Y ¯ P H C A − 1 Y ¯ p , = L B + L γ i , σ − 2 , B , where one frame consists of N consecutive OFDM blocks and the channel length is M . L B ∝ M l o g B + T r B − 1 Y ¯ P H C A i + 1 − 1 Y ¯ p , L γ i , σ − 2 , B ∝ N l o g C A + T r B i − 1 Y ¯ P H C A − 1 Y ¯ p . (4)

Through the optimization of L , we can obtain the updated formula of parameter estimation. First, we calculate the partial derivative of L B to update the formula of the time domain correlation matrix B : ∂ L B ∂ B = M B − 1 − B − 1 Y ¯ P H C A − 1 B − 1 .

We adopt the learning rule for B using ∂ L B ∂ B = 0 : B = Y ¯ P H C A − 1 Y ¯ p M .

Furthermore, the noise variance σ 2 is obtained using the null subcarriers. We rewrite Equation 4 using the Woodbury formula to obtain the model as follows: L γ i = N l o g C − i + T r B − 1 Y ¯ P H C − i − 1 Y ¯ p + N l o g 1 + γ i s i − T r B − 1 γ − 1 + s i q i H q i , = L − i + L i , where L − i ≜ N l o g C − i + T r B − 1 Y ¯ P H C − i − 1 Y ¯ p , (5) L i ≜ log 1 + γ i s i − T r B − 1 γ − 1 + s i q i H q i , (6)

which only depends on γ i . From Equations 5, 6, we can obtain the following: C − i ≜ σ 2 I + ∑ j , j ≠ i M γ j Φ j Φ j H , s i ≜ Φ i H C − i − 1 Φ i , q i ≜ Φ i H C − i − 1 Y .

Setting ∂ L i ∂ γ i = 0 , we have the following updated rule: γ i = q i B − 1 q i H / N − s i s i 2 . (7)

It should be noted that the updated result obtained from Equation 7 is a scalar, where q i B − 1 q i H is q i B 2 , which is the norm of the covariance matrix B .

The steps of the FTMSBL algorithm are provided in Table 2.

The computational complexity of the TMSBL algorithm mainly consists of covariance Σ and mean. The covariance can be solved with marginal cost O M 3 .The computational complexity of mean is O N P M 2 .

The computational complexity of the proposed method is much lower than that of FMLM optimization and Woodbury decomposition. The main calculation quantities are covariance Σ and mean μ . Σ = σ − 2 Φ P H Φ P + Γ r − 1 − 1 , = Γ − Γ Φ P H C A − 1 Φ P Γ , where the covariance Σ can be solved using cost O P 3 + P M 2 , and the total cost of formula μ is O N P M 2 . Therefore, the complexity of the FTMSBL algorithm for one iteration is O P 3 + N P M 2 . The average of the statistical results of algorithm time consumption were calculated and are provided in Table 3.

Table 3 shows the average time consumption of different algorithms. All algorithms are simulated using an Intel (R) Core (TM) i7-8750H CPU@2.2GHz. In terms of computation time, the LS algorithm has very high computational efficiency and the shortest time consumption among all algorithms. The average simulation time of the OMP algorithm is 10 times more than that of LS. Among the SBL algorithms, the FTMSBL algorithm takes the shortest time, approximately 12 s, while the rest of the algorithms take far more than 20 s. This is mainly because the SBL algorithm has a relatively high number of iterations under high accuracy requirements, followed by the high complexity of calculating complex measurement matrices and solving the inverse of matrices during each iteration process. The FTMSBL algorithm essentially avoids the partial matrix inversion process in the traditional SBL algorithm, effectively shortening the computation time.