CN103095639A - Ofdm水声通信并行迭代ici消除方法 - Google Patents

Abstract

本发明的目的在于提供OFDM水声通信并行迭代ICI消除方法，包括以下步骤：建立OFDM水声通信系统，建立水声时变多径信道模型，测量系统中的多普勒因子，估计时变信道的频率响应，估计ICI分量及信号干扰噪声比SINR，对接收信号进行并行迭代MMSE均衡。本发明将线性近似法与频域最小均方误差法结合起来，通过迭代改善性能，从而提出一种并行迭代ICI消除方法并将其引入OFDM水声通信系统中，能有效对抗水声时变信道下产生的载波间干扰，对信道时变速度变化反应稳健。避免了高阶矩阵的求逆运算，降低了方法复杂度，提高了运算速度。

Description

OFDM水声通信并行迭代ICI消除方法

技术领域

本发明涉及的是一种水声移动通信的载波间干扰抑制方法。

背景技术

在复杂多变的水声信道中，尤其在移动通信的条件下，多普勒效应连同多径时延引起的时频选择性衰落，定时，载波和采样频偏等因素，都会打破OFDM系统中各路子载波间的无干扰并行传输，产生ICI，影响数据传输质量。其中尤其以多普勒拓展引起的ICI太复杂而难于进行消除。因而在水声移动通信中，子载波干扰抑制技术成为了决定系统性能的关键技术，日益引起水声领域国内外研究者的广泛关注。

在无线电通信领域目前已有很多行之有效的ICI消除方法。但这些方法大多只针对由频率偏移引起的ICI。在时变信道下，面对由多普勒拓展引起的ICI时常常会束手无策。由于声波在海水中的传播速度与无线电中光速相差好几个量级，产生的多普勒拓展比无线电中大得多。而且相对来说，水声通信系统工作在低频段，带宽也要窄得多，OFDM子载波间隔很小。因而，水声通信中不能简单地将多普勒频偏当成一致频偏来处理，无线电通信领域中的大多数多普勒补偿和ICI抑制方法并不能直接移植到水声通信中。

水声通信中对子载波干扰抑制技术的研究并不多。尤其是时变信道下的ICI抑制技术还很不成熟，通常被当成多普勒问题来看待，并没有特别的针对其进行处理的方法。一般地，在OFDM水声通信中，频偏引起的ICI的消除被归为同步问题，从同步的角度考虑如何有效消除。信道时变引起的ICI则往往被当成多普勒问题进行处理。多普勒估计大多基于模糊函数，补偿则多利用变采样原理。多普勒估计主要有两类，多谱勒频移预估计及块估计方法。常用的包括，插入CW脉冲或者单频信号作为测频序列，在OFDM帧首尾加入同步信号（例如LFM）等。这些估计方法对于固定的多普勒频偏能有效估计，但对于非一致多普勒估计精度却有限。估计出多普勒后可通过插值处理，利用变采样滤波器或高分辨率DFT方法进行补偿。通常利用上述方法对多普勒进行初步补偿后，还存在残留的ICI。为进一步改善性能，事实上还必须对残余的ICI进行消除。而目前在OFDM水声通信领域中，针对两步多普勒补偿的研究还很少。

在水声移动通信中，要消除由信道时变引起的ICI，如果不从多普勒补偿的角度，最好的方式就是进行时变信道估计，然后在接收端进行ICI消除后再均衡。而关系到ICI消除效果的关键步骤就是时变信道的估计。如果借助一个OFDM符号时间内的线性近似模型来进行时变信道估计，仅需分别估计出包含信道时不变及时变信息的两组向量就能获取完整的信道冲击响应。

在水声领域目前还未有相关的研究。在无线电通信中，时变信道线性近似方面的研究已有一些。文献（M.Johnson,L.Freitag and M.Stojanovic.ImprovedDoppler tracking and correction for underwater acoustic communications[C]//in Proc.ICASSP’97.Munich,Germany，1997:575–578.）就利用了上述思想，基于最小平方（least square，LS）准则获取信道时变信息后通过多项式近似法估计出ICI并将其消除，再通过迭代进一步提高性能。文献（Bayan S.Sharif，Jeff Neasham,Oliver R.Hinton and AlanE.Adams.A Computationally Efficient DopplerCompensation System for Underwater Acoustic Communications[J].OceanicEngineering,2000,25(1):52-61.）同样基于分段线性近似模型，采用期望最大化(estimation maximum，EM)迭代方法来提高符号平均信道脉冲响应的估计精度，并在迭代过程中进行带状子载波间干扰抑制。文献（Sean F.Mason,Christian R.Berger,Shengli Zhou and Peter Willett.Detection,Synchronization,and DopplerScale Estimation with Multicarrier Waveforms in Underwater AcousticCommunication[J].Journal on selected areas in communications,2008,26(9):1638-1649.）先使用MMSE方法对传输信号进行初值估计，再通过迭代干扰抵消和均衡相结合，来消除ICI。文献（Maja Sliskovic.Sampling Frequency OffsetEstimation and Correction in OFDM Systems[C]//in Proc.ICECS 2001.Malta,2001:437-440.）则利用信道译码输出软信息来当作信息比特的先验信息，再根据子载波间干扰矩阵来消除干扰分量，经过若干次迭代来达到性能优化。这两种方法利用了迭代ICI消除与均衡相结合的思路，但在ICI消除环节，使用的却是更新均衡因子或干扰自消除的办法。

发明内容

本发明的目的在于提供能有效对抗水声时变信道下产生的载波间干扰的OFDM水声通信并行迭代ICI消除方法。

本发明的目的是这样实现的：

本发明OFDM水声通信并行迭代ICI消除方法，其特征是：

（1）建立OFDM水声通信系统；

（2）建立水声时变多径信道模型；

（3）测量系统中的多普勒因子；

（4）估计时变信道的频率响应；

（5）估计ICI分量及信号干扰噪声比SINR；

（6）对接收信号进行并行迭代MMSE均衡。

本发明还可以包括：

1、用WSSUS模型来描述水声时变多径信道，冲击响应的相位服从均匀分布，能量或者包络呈Rayleigh分布，冲激响应用抽头延迟线模型来实现，各路径的时延固定不变，冲击响应h(n，λ)的表达式为：

$h(n,\mathrm{\λ})=\underset{l=0}{\overset{L-1}{\mathrm{\Σ}}}{h}_l\left(n\right)\mathrm{\δ}(\mathrm{\λ}-l)=\underset{l=0}{\overset{L-1}{\mathrm{\Σ}}}{h}_l{e}^{j2\mathrm{\π}{\mathrm{\ϵ}}_{l}n/N}\·\mathrm{\δ}(\mathrm{\λ}-l),0\≤n\≤N-1$

其中N表示总子载波数，n为子载波序号，λ表示时延长度，ε_l＝f_dlT_s是归一化多普勒频偏，δ(*)为冲击函数，h_l(n)、f_dl分别是第l个抽头上的加权系数和多普勒频移，L为抽头延迟线的长度，各抽头间隔相等，均为采样时间T_s，WSSUS模型采用Jakes模型来实现，所有信道抽头都是独立的，并且每个抽头都是一个独立的均值为零的广义平稳高斯随机过程，其信道的时间差相关函数R_h(Δn)为：

$R_h\left(\mathrm{\Δn}\right)=E\{h_l\left(n\right)*h_l^{*}(n+\mathrm{\Δn})\}={\mathrm{\δ}}_l^2{J}_0\left({2\mathrm{\πf}}_{d\max }\right|\mathrm{\Δn}\left|T_s\right)$

E{*}表示取均值，δ_l为对应各路径时延的冲击函数，Δn为子载波序号差，J₀(*)为零阶第一类Bessel函数，f_dmax是最大多普勒频偏。

2、采用多普勒块估计方法对多普勒因子进行测量，向发射OFDM帧首尾加入同步信号LFM，在接收端通过同步头检测获取每帧OFDM信号的LEM相关峰，定义Peak_i为第i帧OFDM信号的LEM相关峰位置，下标i表示OFDM信号帧序号，根据原信号帧长度L_SignalFrame，获取每帧信号上的多普勒因子 ${\mathrm{\α}}_i=\frac{{\mathrm{Peak}}_{i+1}-{\mathrm{Peak}}_i}{L_{\mathrm{SignalFrame}}}-1.$

3、时变信道频率响应的估计基于分段线性近似信道模型，并通过梳状导频辅助的信道估计方式来实现：

（1）估计中间时刻的信道冲击响应

定义新参数

表示第l条路径上的信道冲击响应在一个OFDM符号时间内的时间平均值，令υ＝ceil(τ_max,T_s)表示最大归一化时延，则υ+1≤G，G为循环前缀数目，τ_max为信道最大多径时延，然后等间距地向OFDM符号中插入梳状导频，导频数用M表示，用H(m_ii，m_ii)表示整个OFDM符号中m_ii位置处的导频序列上的信道频率响应，m_ii＝ii×N/M，0≤ii≤M-1，在设置系统参数时令M等于G的整数倍，并且确保M的值可表示为2^k，得到信道冲击响应平均值

的表达式：

$h_l^{\mathrm{ave}}=\frac{1}{M}\underset{\mathrm{ii}=0}{\overset{M-1}{\mathrm{\Σ}}}H(m_{\mathrm{ii}},m_{\mathrm{ii}})\·e^{j2\mathrm{\πii\μl}/M}$

当n＝N/2-1时，

取得最小值，即在OFDM符号中间时刻上的信道冲击响应

近似等于整个OFDM符号时间内信道冲击响应时间平均值即：

$h_l=(\frac{N}{2}-1)\≈h_l^{\mathrm{ave}};$

（2）估计信道冲击响应时变斜率

以OFDM符号中间时刻为基准，将一个OFDM符号划为两个区域，每个OFDM符号在区域1和区域2上的各存在一组信道变化斜率向量，分别用α^r1和α^r2表示，即 ${\mathrm{\α}}^{\mathrm{rq}}={(a_0^{\mathrm{rq}},a_1^{\mathrm{rq}},...,a_{L-1}^{\mathrm{rq}})}^T,q=1\mathrm{or}2,$ 其中q表示区域，下标表示路径数，令u为当前OFDM符号序号，T、T_g分别表示OFDM符号长度和循环前缀长度，则：

${\mathrm{\α}}_l^{\mathrm{rq}}=\frac{h_l^{u+q-1}(\frac{N}{2}-1)-h_l^{u+q-2}(\frac{N}{2}-1)}{T+T_g},0\≤l\≤L-1,q=1\mathrm{or}2$

利用

α^r1和α^r2进行一阶线性插值逼近，获取当前OFDM符号上的完整信道冲击响应

$h_l^u\left(n\right)=\left\{\begin{array}{c}{h}_l(\frac{N}{2}-1)+(n+1-\frac{N}{2})\×{\mathrm{\α}}_l^{r1}\×T_s,0\≤n\≤\frac{N}{2}-1\\ h_l(\frac{N}{2}-1)+(n+1-\frac{N}{2})\×{\mathrm{\α}}_l^{r2}\×T_s,\frac{N}{2}\≤n\≤N-1\end{array}\right.;$

（3）获取信道频率响应矩阵

信道频率响应H表示为：

$H=H_{\mathrm{mid}}+C^{r1}\×H_{\mathrm{slope}}^{r1}+C^{r2}\×H_{\mathrm{slope}}^{r2}---\left(7\right)$

其中， $H_{\mathrm{mid}}=\mathrm{diag}\left\{\mathrm{FFT}\left(\right[h_0(\frac{N}{2}-1)h_1(\frac{N}{2}-1)...h_{G-1}(\frac{N}{2}-1)0...0\left]\right)\right\},$

$H_{\mathrm{slope}}^{\mathrm{rq}}=\mathrm{diag}\left\{\mathrm{FFT}\left(\left[\begin{array}{cccc}{\mathrm{\α}}_0^{\mathrm{rq}}& {\mathrm{\α}}_1^{\mathrm{rq}}& \·\·\·{\mathrm{\α}}_{G-1}^{\mathrm{rq}}& 0\·\·\·0\end{array}\right]\right)\right\},$

$C^{r1}(k,m)=\mathrm{Ts}\×\left\{\begin{array}{cc}\frac{-0.5}{{1-e}^{\frac{-j2\mathrm{\π}(k-m)}{N}}}+\frac{1-{(-1)}^{k-m}}{{(1-e^{\frac{-j2\mathrm{\π}(k-m)}{N}})}^2\×N}& k\≠m\\ -\frac{1}{4}-\frac{N}{8}& k=m\end{array}\right.,$

$C^{r2}(k,m)=\mathrm{Ts}\×\left\{\begin{array}{cc}\frac{-0.5}{1-e^{\frac{-j2\mathrm{\π}(k-m)}{N}}}-\frac{1-{(-1)}^{k-m}}{{(1-e^{\frac{-j2\mathrm{\π}(k-m)}{N}})}^2\×N}& k\≠m\\ -\frac{1}{4}+\frac{N}{8}& k=m\end{array}\right.,0\≤k,m\≤N-1.$

4、所述的ICI分量及信号干扰噪声比SINR的估计方法为：

在接收端去除CP后的信号y(n)可表示为：

$y\left(n\right)=s\left(n\right)\⊗h(n,\mathrm{\λ})+w\left(n\right),0\≤n\≤N-1$

其中s(n)为发送信号，w(n)为噪声序列，将y(n)进行DFT变换，得到频域接收信号Y(k)的表达式，其中S(k)和W(k)分别表示信号及噪声的频域表达式，ICI(k)为ICI的表达式：

$Y\left(k\right)=H\left(k\right)S\left(k\right)+\underset{m\≠k}{\overset{N-1}{\underset{m=0}{\mathrm{\Σ}}}}H(k,m)S\left(m\right)+W\left(k\right)$

$=H(k,k)S\left(k\right)+\mathrm{ICI}\left(k\right)+W\left(k\right)$

再结合前一次均衡所估计出的发送信号p表示迭代次数，完成对ICI的估计：

$\mathrm{ICI}=H^{\′\′}{\hat{S}}^{(p-1)}$

其中H'＝diag(H)，H″＝H-H'；

结合ICI表达式及Jakes模型下信道的时间差相关函数获得信号干扰噪声比SINR：

$\mathrm{SINR}={E\left|S\left(m\right)\right|}^2\frac{\frac{1}{N^2}[N+2\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}(N-n)J_0\left({2\mathrm{\πf}}_{d\max }{T}_{s}n\right)]}{1-\frac{1}{N^2}[N+2\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}(N-n)J_0\left({2\mathrm{\πf}}_{d\max }{T}_{s}n\right)]+\frac{{E\left|W\left(m\right)\right|}^2}{{E\left|S\left(m\right)\right|}^2}}$

其中E|S(m)|²、E|W(m)|²分别表示信号及噪声的功率，最大多普勒频偏f_dmax通过f_dmax＝αf_H得到，α为多普勒因子,f_H为信号频率上限。

5、所述的并行迭代MMSE均衡为：

（1）对MMSE准则进行修正

将均衡矩阵H_mmse中的信噪比SNR用信号干扰噪声比SINR代替，即

$H_{\mathrm{mmse}}=R_{H_{\mathrm{ls}}{H}_{\mathrm{ls}}}{(R_{H_{\mathrm{ls}}{H}_{\mathrm{ls}}}+\frac{\mathrm{\β}}{\mathrm{SINR}}I)}^{-1}{H}_{\mathrm{ls}}$

其中H_ls是LS估计得到的信道频响，

是其自相关函数，I是单位矩阵，对于QPSK调制β取1；

（2）进行并行迭代MMSE均衡

ξ为迭代终止阈值，p_max表示最大迭代次数；

①求出在迭代均衡过程中固定不变的矩阵C^r1，C^r2及信号干扰噪声比SINR；

②进行LS估计，获得导频位置处的信道频响并进行一阶线性频域插值获得所有子载波上的信道频响H_ls，然后借助步骤①中求得的SINR，得到H_mmse；再由

完成MMSE均衡，得到发送信号的估计值

③初始化迭代次数p＝1，令第一次迭代前的接收信号值Y⁽⁰⁾＝Y，Y为接收符号，第一次迭代前的接收信号值估计值

④估计出导频位置处的信道频响

$\hat{H}(m_{\mathrm{ii}},m_{\mathrm{ii}})=\frac{Y^{p-1}\left(m_{\mathrm{ii}}\right)}{P\left(m_{\mathrm{ii}}\right)}=H(m_{\mathrm{ii}},m_{\mathrm{ii}})+\frac{\mathrm{ICI}\left(m_{\mathrm{ii}}\right)+W\left(m_{\mathrm{ii}}\right)}{P\left(m_{\mathrm{ii}}\right)}$

其中P表示导频序列；

⑤获取中间时刻信道冲击响应h^ave，继而得出信道时变斜率α^r1，α^r2；

⑥获取矩阵H_mid，

及

⑦求出矩阵H，H'及H″，并结合前一次估计出的

估计出接收信号中的ICI^(p)分量；

⑧从接收符号Y中将ICI分量减掉，令Y^p＝Y-ICI^(p)；

⑨获得H_mmse再由

完成MMSE均衡，得到p次迭代后的发送信号的估计值

⑩当估计值足够收敛即满足

或p≥p_max时，终止迭代；否则迭代次数p增1，返回到步骤④。

6、p_max小于或等于8。

本发明的优势在于：本发明将线性近似法与频域最小均方误差法结合起来，通过迭代改善性能，从而提出一种并行迭代ICI消除方法并将其引入OFDM水声通信系统中，能有效对抗水声时变信道下产生的载波间干扰，对信道时变速度变化反应稳健。避免了高阶矩阵的求逆运算，降低了方法复杂度，提高了运算速度。

附图说明

图1为本发明框图；

图2为本发明的流程图；

图3为本发明水声通信系统原理框图；

图4为本发明近似线性信道模型的示意图；

图5为时变信道估计及ICI估计部分流程图。

具体实施方式

下面结合附图举例对本发明做更详细地描述：

结合图1～5，基于信道在一个OFDM符号时间内近似线性变化的假设，在不改变OFDM系统结构的情况下，利用梳状导频辅助估计出每个OFDM符号周期内中间时刻的信道冲击响应，同时从相邻OFDM符号上的信道信息中提取出信道冲击响应的时变斜率。然后，借助得到的这两组向量进行分段线性插值，进一步估计出其他时刻的信道冲击响应，最终获取整个OFDM符号时间内的信道时变信息。接下来利用推导得出的ICI表达式，结合前一次均衡得到的发送符号的估计值，就可以得到接收信号中ICI的估计值，进而将其从中减掉。在此基础上再次进行均衡，得到发送符号的新的估计值。

本发明分为以下步骤：

1、建立OFDM水声通信系统

图3给出了一个典型的OFDM水声通信系统原理框图。OFDM水声通信系统离不开图中几个基本环节。考虑到水声信道的多途效应和随机时变，多普勒影响，水下噪声，实验环境，硬件性能等因素，图像传输时往往会增加交织等步骤，通常还不得不考虑同步等问题，可能还需进行多普勒补偿。

2、建立水声时变多径信道模型

水声信道是时变多径信道，不同传输路径上的时延和衰减系数都是时变的。对于这样复杂的信道，通常可以用WSSUS模型来描述。在这个模型下，冲击响应的相位服从均匀分布，能量或者包络呈Rayleigh分布。离散WSSUS信道的冲激响应往往可以用抽头延迟线模型来实现。此时我们通常认为各路径的时延是固定不变的，从而减少信道建模的复杂度。其冲击响应h(n,λ)的表达式为：

$h(n,\mathrm{\λ})=\underset{l=0}{\overset{L-1}{\mathrm{\Σ}}}{h}_l\left(n\right)\mathrm{\δ}(\mathrm{\λ}-l)=\underset{l=0}{\overset{L=1}{\mathrm{\Σ}}}{h}_l{e}^{j2\mathrm{\π}{\mathrm{\ϵ}}_{l}n/N}\·\mathrm{\δ}(\mathrm{\λ}-l)---\left(1\right)$

其中ε_l＝f_dlT_s是归一化多普勒频偏，δ(*)为冲击函数。h_l(n)，f_dl分别是第l个抽头上的加权系数和多普勒频移。f_dl的存在导致了h_l(n)的时变性。抽头延迟线的长度L由信道最大多径时延τ_max决定，各抽头间隔相等，均为采样时间T_s。在Jakes模型中，所有信道抽头都是独立的，并且每个抽头都是一个独立的均值为零的广义平稳高斯随机过程。其信道的时间差相关函数R_h(Δn)为：

$R_h\left(\mathrm{\Δn}\right)=E\{h_l\left(n\right)*h_l^{*}(n+\mathrm{\Δn})\}={\mathrm{\δ}}_l^2{J}_0\left({2\mathrm{\πf}}_{d\max }\right|\mathrm{\Δn}\left|T_s\right)---\left(2\right)$

在上式中，E{*}表示取均值，δ_l为对应各路径时延的冲击函数。J₀(*)为零阶第一类Bessel函数，f_dmax是最大多普勒频偏。

3、测量系统中的多普勒因子

多普勒因子的估计方法主要有两类，多谱勒频移预估计及块估计方法。为降低方法的复杂度，我们采用多普勒块估计方法对多普勒因子进行测量。首先，向发射OFDM帧首尾加入同步信号LFM。在接收端通过同步头检测可获取每帧OFDM信号的LEM相关峰。最后借助相邻相关峰间距Peak_i+1-Peak_i及原信号帧长度L_SignalFrame即可获取每帧信号上的多普勒因子

4、估计时变信道的频率响应

本方法涉及的时变信道估计主要基于分段线性近似信道模型，并通过梳状导频辅助的信道估计方式来实现。

（1）估计中间时刻的信道冲击响应

首先定义一个新的参数用来表示第l条路径上的信道冲击响应在一个OFDM符号时间内的时间平均值。下面考虑梳状导频辅助下的时变信道估计。

令υ＝ceil(τ_max,Ts)，表示最大归一化时延。那么，υ+1≤G，其中G为循环前缀数目。然后，我们等间距地向OFDM符号中插入梳状导频，导频数用M表示。用H(m_ii，m_ii)表示整个OFDM符号中m_ii位置处的导频序列上的信道频率响应，其中，m_ii=ii×N/M，0≤ii≤M-1。经过一系列公式推导发现H(m_ii，m_ii)其实就是信道冲击响应平均值

的傅里叶变换，相当于对长为G的

进行补零操作后再进行M点的DFT。此时要对

进行估计，则只需进行IDFT操作。而根据DFT插值的原理，要想通过H(m_ii，m_ii)准确无误地恢复出

通常需要在设置系统参数时令M等于G的整数倍（假定比值为μ），并且确保M的值可表示为2^k。最后得到信道冲击响应平均值

的表达式：

$h_l^{\mathrm{ave}}=\frac{1}{M}\underset{\mathrm{ii}=0}{\overset{M-1}{\mathrm{\Σ}}}H(m_{\mathrm{ii}},m_{\mathrm{ii}})\·e^{j2\mathrm{\πii\μl}/M}---\left(3\right)$

通过讨论知道对于第l条路径而言，当n＝N/2-1时，

取得最小值，即在OFDM符号中间时刻上的信道冲击响应近似等于整个OFDM符号时间内信道冲击响应时间平均值

如式(4)所示，这也与线性信道模型的假设相吻合。

$h_l(\frac{N}{2}-1)\≈h_l^{\mathrm{ave}}---\left(4\right)$

（2）估计信道冲击响应时变斜率

基于线性信道模型^[8]，通常我们认为信道在一个OFDM符号时间内的变化斜率是固定不变的。本方法以OFDM符号中间时刻为基准，将一个OFDM符号划为两个区域来探究信道冲击响应时变斜率的估计方法。图4给出了所建立的近似线性信道模型的示意图。

由上图可以知道每个OFDM符号在区域1和区域2上的各存在一组信道变化斜率向量，分别用α^r1和α^r2表示，即其中q表示区域，下标表示路径数。下面对α^r1和α^r2进行求解。T,T_g分别表示OFDM符号长度和循环前缀长度。

${\mathrm{\α}}_l^{\mathrm{rq}}=\frac{h_l^{i+q-1}(\frac{N}{2}-1)-h_l^{i+q-2}(\frac{N}{2}-1)}{T+T_g},0\≤l\≤L-1,q=1\mathrm{or}2---\left(5\right)$

一旦中间时刻的信道冲击响应被准确无误地估计出来，就很容易求出两组信道变化斜率向量。接下来要做的就是利用

α^r1和α^r2进行一阶线性插值逼近，获取当前OFDM符号上的完整信道冲击响应。

$h_l^i\left(n\right)=\left\{\begin{array}{c}{h}_l(\frac{N}{2}-1)+(n+1-\frac{N}{2})\×{\mathrm{\α}}_l^{r1}\×T_s,0\≤n\≤\frac{N}{1}-1\\ h_l(\frac{N}{2}-1)+(n+1-\frac{N}{2})\×{\mathrm{\α}}_l^{r2}\×T_s,\frac{N}{2}\≤n\≤N-1\end{array}\right.---\left(6\right)$

（3）获取信道频率响应矩阵

事实上，只要推导得出时变信道的时域冲击响应函数，进行FFT变换到频域就能直接获取信道频率响应。考虑到插值及相应FFT变换的算法复杂度，本发明经过公式推导得到了信道频率响应的闭式表达式，从而仅由中间时刻的信道冲击响应及信道变化斜率向量实现了信道频率响应矩阵的准确估计。

仍然针对图4所示的近似线性信道模型进行讨论，主要基于时域信道矩阵h的循环移位Toeplitz特性进行推导，最终可将信道频率响应H表示为：

$H=H_{\mathrm{mid}}+C^{r1}\×H_{\mathrm{slope}}^{r1}+C^{r2}\×H_{\mathrm{slope}}^{r2}---\left(7\right)$

其中，

$H_{\mathrm{mid}}=\mathrm{diag}\left\{\mathrm{FFT}\left(\right[h_0(\frac{N}{2}-1h_1(\frac{N}{2}-1)...h_{G-1}(\frac{N}{2}-1)0...0])\right\}---\left(8\right)$

$H_{\mathrm{slope}}^{\mathrm{rq}}=\mathrm{diag}\left\{\mathrm{FFT}\left(\left[\begin{array}{cccc}{\mathrm{\α}}_0^{\mathrm{rq}}& {\mathrm{\α}}_1^{\mathrm{rq}}& \·\·\·{\mathrm{\α}}_{G-1}^{\mathrm{rq}}& 0\·\·\·0\end{array}\right]\right)\right\}---\left(9\right)$

C^r1，C^r2由时域信道矩阵h的循环移位Toeplitz特性导出，推导过程较为复杂。定义1≤k,m≤N，则其表达式如下：

$C^{r1}(k,m)=\mathrm{Ts}\×\left\{\begin{array}{cc}\frac{-0.5}{1-e^{\frac{-j2\mathrm{\π}(k-m)}{N}}}+\frac{1-{(-1)}^{k-m}}{{(1-e^{\frac{-j2\mathrm{\π}(k-n)}{N}})}^2\×N}& k\≠m\\ -\frac{1}{4}-\frac{N}{8}& k=m\end{array}\right.---\left(10\right)$

$C^{r2}(k,m)=\mathrm{Ts}\×\left\{\begin{array}{cc}\frac{-0.5}{1-e^{\frac{-j2\mathrm{\π}(k-m)}{N}}}-\frac{1-{(-1)}^{k-m}}{{(1-e^{\frac{-j2\mathrm{\π}(k-m)}{N}})}^2\×N}& k\≠m\\ -\frac{1}{4}+\frac{N}{8}& k=m\end{array}\right.---\left(11\right)$

5、估计ICI分量及信号干扰噪声比SINR

时变信道下，传输信号经过时频双选择性信道，同时还存在加性高斯白噪声。于是，在接收端去除CP后的信号y(n)可表示为：

$y\left(n\right)=s\left(n\right)\⊗h(n,\mathrm{\λ})+w\left(n\right),0\≤n\≤N-1---\left(12\right)$

其中，s(n)为发送信号，w(n)为噪声序列。结合离散WSSUS信道的冲激响应表达式，并将y(n)进行DFT变换，可得到频域接收信号Y(k)的表达式。其中，S(k)和W(k)分别表示信号及噪声的频域表达式。为满足ICI分析的需要，可以将其表示为如下形式：

此时再结合前一次均衡所估计出的发送信号就能完成对ICI的估计。

$\mathrm{ICI}=H^{\′\′}{\hat{S}}^{(p-1)}---\left(14\right)$

其中，H'＝diag(H)，H″＝H-H'。接下来只要结合ICI表达式及Jakes模型下信道的时间差相关函数就能推导出信号干扰噪声比SINR。

$\mathrm{SINR}={E\left|S\left(m\right)\right|}^2\frac{\frac{1}{N^2}[N+2\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}(N-n)J_0\left({2\mathrm{\πf}}_{d\max }{T}_{s}n\right)]}{1-\frac{1}{N^2}[N+2\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}(N-n)J_0\left({2\mathrm{\πf}}_{d\max }{T}_{s}n\right)]+\frac{{E\left|W\left(m\right)\right|}^2}{{E\left|S\left(m\right)\right|}^2}}---\left(15\right)$

其中，E|S(m)|²，E|W(m)|²分别表示信号及噪声的功率。最大多普勒频偏f_dmax可由f_dmax＝αf_H计算得到。

6、对接收信号进行并行迭代MMSE均衡

（1）对MMSE准则进行修正

在时变信道下，由于ICI的影响若要确保MMSE均衡的有效性必须对原有的MMSE准则进行修正。在近似线性信道模型下，ICI可以看成是一种加性干扰。在对其进行准确估计后，可以直接从接收信号中将其减去。因此，在基于MMSE准则的并行迭代均衡方法中，可将均衡矩阵中的信噪比SNR用信号干扰噪声比SINR代替，即

$H_{\mathrm{mmse}}=R_{H_{\mathrm{ls}}{H}_{\mathrm{ls}}}{(R_{H_{\mathrm{ls}}{H}_{\mathrm{ls}}}+\frac{\mathrm{\β}}{\mathrm{SINR}}I)}^{-1}{H}_{\mathrm{ls}}---\left(16\right)$

其中，H_ls是LS估计得到的信道频响，

是其自相关函数。I是单位矩阵。对于QPSK调制β取1。

（2）进行并行迭代MMSE均衡

ξ为迭代终止阈值，是预先设置好的。p_max表示最大迭代次数，由于多次迭代势必占据大量时间，一般不宜取得过大。图5给出了时变信道估计及ICI估计部分流程图。

步骤一：求出在迭代均衡过程中固定不变的矩阵C^r1，C^r2及信号干扰噪声比SINR。

步骤二：进行LS估计，获得导频位置处的信道频响

并进行一阶线性频域插值获得所有子载波上的信道频响H_ls。然后借助步骤一中求得的SINR，利用式(16)，得到H_mmse；再由

完成MMSE均衡，得到发送信号的估计值

步骤三：初始化迭代次数p＝1，令第一次迭代前的接收信号值Y⁽⁰⁾＝Y，第一次迭代前的接收信号值估计值

步骤四：估计出导频位置处的信道频响

为下一步求解中间时刻信道冲击响应向量及时变斜率向量做准备，

其中，P表示导频序列。实际上是进行了一次LS估计。其后通过一阶线性频域插值得到整个OFDM符号上的信道频响估计值H_ls。

步骤五：利用式(3-4)，获取中间时刻信道冲击响应h^ave，继而得出信道时变斜率α^r1，α^r2。

步骤六：利用式(8)及(9)，获取矩阵H_mid，

及

步骤七：求出矩阵H，H'及H″，并结合前一次估计出的

估计出接收信号中的ICI^(p)分量。

步骤八：从接收符号Y中将ICI分量减掉，令Y^p＝Y-ICI^(p)。

步骤九：利用式(16)，得到H_mmse再由

完成MMSE均衡，得到发送信号的估计值

步骤十：迭代终止条件判决。当估计值足够收敛即满足

或p≥p_max时，终止迭代，否则迭代次数p增1。返回到步骤四。

如此进行反复迭代，可以不断改善ICI的估计精度，在此基础上进行ICI消除并再次均衡可以获得更好的性能。通常迭代次数越高，其精度越高。但也存在一定的收敛性，一般迭代7-8次，将达到方法性能极限。本方法中通过预先设置一个迭代终止阈值ξ，来控制迭代的终止。再结合具体要求合理设置迭代次数上限p_max，在精度和速度上做一下平衡，从而使本方法达到其最佳性能。

Claims (7)

1.OFDM水声通信并行迭代ICI消除方法，其特征是：

（1）建立OFDM水声通信系统；

（2）建立水声时变多径信道模型；

（3）测量系统中的多普勒因子；

（4）估计时变信道的频率响应；

（5）估计ICI分量及信号干扰噪声比SINR；

（6）对接收信号进行并行迭代MMSE均衡。

2.根据权利要求1所述的OFDM水声通信并行迭代ICI消除方法，其特征是：用WSSUS模型来描述水声时变多径信道，冲击响应的相位服从均匀分布，能量或者包络呈Rayleigh分布，冲激响应用抽头延迟线模型来实现，各路径的时延固定不变，冲击响应h(n,λ)的表达式为：

$h(n,\mathrm{\λ})=\underset{l=0}{\overset{L-1}{\mathrm{\Σ}}}{h}_l\left(n\right)\mathrm{\δ}(\mathrm{\λ}-l)=\underset{l=0}{\overset{L-1}{\mathrm{\Σ}}}{h}_l{e}^{j2\mathrm{\π}{\mathrm{\ϵ}}_{l}n/N}\·\mathrm{\δ}(\mathrm{\λ}-l),0\≤n\≤N-1$

其中N表示总子载波数，n为子载波序号，λ表示时延长度，ε_l＝f_dlT_s是归一化多普勒频偏，δ(*)为冲击函数，h_l(n)、f_dl分别是第l个抽头上的加权系数和多普勒频移，L为抽头延迟线的长度，各抽头间隔相等，均为采样时间T_s，WSSUS模型采用Jakes模型来实现，所有信道抽头都是独立的，并且每个抽头都是一个独立的均值为零的广义平稳高斯随机过程，其信道的时间差相关函数R_h(Δn)为：

$R_h\left(\mathrm{\Δn}\right)=E\{h_l\left(n\right)*h_l^{*}(n+\mathrm{\Δn})\}={\mathrm{\δ}}_l^2{J}_0\left({2\mathrm{\πf}}_{d\max }\right|\mathrm{\Δn}\left|T_s\right)$

E{*}表示取均值，δ_l为对应各路径时延的冲击函数，Δn为子载波序号差，J₀(*)为零阶第一类Bessel函数，f_dmax是最大多普勒频偏。

3.根据权利要求1所述的OFDM水声通信并行迭代ICI消除方法，其特征是：采用多普勒块估计方法对多普勒因子进行测量，向发射OFDM帧首尾加入同步信号LFM，在接收端通过同步头检测获取每帧OFDM信号的LEM相关峰，定义Peak_i为第i帧OFDM信号的LEM相关峰位置，下标i表示OFDM信号帧序号，根据原信号帧长度L_SignalFrame，获取每帧信号上的多普勒因子 ${\mathrm{\α}}_i=\frac{{\mathrm{Peak}}_{i+1}-{\mathrm{Peak}}_i}{L_{\mathrm{SignalFrame}}}-1.$

4.根据权利要求1所述的OFDM水声通信并行迭代ICI消除方法，其特征是：时变信道频率响应的估计基于分段线性近似信道模型，并通过梳状导频辅助的信道估计方式来实现：

（1）估计中间时刻的信道冲击响应

定义新参数

表示第l条路径上的信道冲击响应在一个OFDM符号时间内的时间平均值，令υ＝ceil(τ_max,T_s)表示最大归一化时延，则υ+1≤G，G为循环前缀数目，τ_max为信道最大多径时延，然后等间距地向OFDM符号中插入梳状导频，导频数用M表示，用H(m_ii，m_ii)表示整个OFDM符号中m_ii位置处的导频序列上的信道频率响应，m_ii＝ii×N/M，0≤ii≤M-1，在设置系统参数时令M等于G的整数倍，并且确保M的值可表示为2^k，得到信道冲击响应平均值

的表达式：

$h_l^{\mathrm{ave}}=\frac{1}{M}\underset{\mathrm{ii}=0}{\overset{M-1}{\mathrm{\Σ}}}H(m_{\mathrm{ii}},m_{\mathrm{ii}})\·e^{j2\mathrm{\πii\μl}/M}$

当n＝N/2-1时，

取得最小值，即在OFDM符号中间时刻上的信道冲击响应

近似等于整个OFDM符号时间内信道冲击响应时间平均值即：

$h_l(\frac{N}{2}-1)\≈h_l^{\mathrm{ave}};$

（2）估计信道冲击响应时变斜率

以OFDM符号中间时刻为基准，将一个OFDM符号划为两个区域，每个OFDM符号在区域1和区域2上的各存在一组信道变化斜率向量，分别用α^r1和α^r2表示，即 ${\mathrm{\α}}^{\mathrm{rq}}={(a_0^{\mathrm{rq}},a_1^{\mathrm{rq}},\·\·\·,a_{L-1}^{\mathrm{rq}})}^T,q=1\mathrm{or}2,$ 其中q表示区域，下标表示路径数，令u为当前OFDM符号序号，T、T_g分别表示OFDM符号长度和循环前缀长度，则：

${\mathrm{\α}}_l^{\mathrm{rq}}=\frac{h_l^{u+q-1}(\frac{N}{2}-1)-h_l^{u+q-2}(\frac{N}{2}-1)}{T+T_g},0\≤l\≤L-1,q=1\mathrm{or}2$

利用

α^r1和α^r2进行一阶线性插值逼近，获取当前OFDM符号上的完整信道冲击响应

$h_l^u\left(n\right)=\left\{\begin{array}{c}{h}_l(\frac{N}{2}-1)+(n+1-\frac{N}{2})\×{\mathrm{\α}}_l^{r1}\×T_s,0\≤n\≤\frac{N}{2}-1\\ h_l(\frac{N}{2}-1)+(n+1-\frac{N}{2})\×{\mathrm{\α}}_l^{r2}\×T_s,\frac{N}{2}\≤n\≤N-1\end{array}\right.;$

（3）获取信道频率响应矩阵

信道频率响应H表示为：

$H=H_{\mathrm{mid}}+C^{r1}\×H_{\mathrm{slope}}^{r1}+C^{r2}\×H_{\mathrm{slope}}^{r2}---\left(7\right)$

其中， $H_{\mathrm{mid}}=\mathrm{diag}\left\{\mathrm{FFT}\left(\right[h_0(\frac{N}{2}-1)h_1(\frac{N}{2}-1)...h_{G-1}(\frac{N}{2}-1)0...0\left]\right)\right\},$

$H_{\mathrm{slope}}^{\mathrm{rq}}=\mathrm{diag}\left\{\mathrm{FFT}\left(\left[\begin{array}{cccc}{\mathrm{\α}}_0^{\mathrm{rq}}& {\mathrm{\α}}_1^{\mathrm{rq}}& \·\·\·{\mathrm{\α}}_{G-1}^{\mathrm{rq}}& 0\·\·\·0\end{array}\right]\right)\right\},$

$C^{r1}(k,m)=\mathrm{Ts}\×\left\{\begin{array}{cc}\frac{-0.5}{{1-e}^{\frac{-j2\mathrm{\π}(k-m)}{N}}}+\frac{1-{(-1)}^{k-m}}{{(1-e^{\frac{-j2\mathrm{\π}(k-m)}{N}})}^2\×N}& k\≠m\\ -\frac{1}{4}-\frac{N}{8}& k=m\end{array}\right.,$

$C^{r2}(k,m)=\mathrm{Ts}\×\left\{\begin{array}{cc}\frac{-0.5}{{1-e}^{\frac{-j2\mathrm{\π}(k-m)}{N}}}-\frac{1-{(-1)}^{k-m}}{{(1-e^{\frac{-j2\mathrm{\π}(k-m)}{N}})}^2\×N}& k\≠m\\ -\frac{1}{4}+\frac{N}{8}& k=m\end{array}\right.,$

$0\≤k,m\≤N-1.$

5.根据权利要求1所述的OFDM水声通信并行迭代ICI消除方法，其特征是：所述的ICI分量及信号干扰噪声比SINR的估计方法为：

在接收端去除CP后的信号y(n)可表示为：

$y\left(n\right)=s\left(n\right)\⊗h(n,\mathrm{\λ})+w\left(n\right),0\≤n\≤N-1$

其中s(n)为发送信号，w(n)为噪声序列，将y(n)进行DFT变换，得到频域接收信号Y(k)的表达式，其中S(k)和W(k)分别表示信号及噪声的频域表达式，ICI(k)为ICI的表达式：

$Y\left(k\right)=H\left(k\right)S\left(k\right)+\underset{m\≠k}{\overset{N-1}{\underset{m=0}{\mathrm{\Σ}}}}H(k,m)S\left(m\right)+W\left(k\right)$

$=H(k,k)S\left(k\right)+\mathrm{ICI}\left(k\right)+W\left(k\right)$

再结合前一次均衡所估计出的发送信号

p表示迭代次数，完成对ICI的估计：

$\mathrm{ICI}=H^{\′\′}{\hat{S}}^{(p-1)}$

其中H'＝diag(H)，H″＝H-H'；

结合ICI表达式及Jakes模型下信道的时间差相关函数获得信号干扰噪声比SINR：

$\mathrm{SINR}={E\left|S\left(m\right)\right|}^2\frac{\frac{1}{N^2}[N+2\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}(N-n)J_0\left({2\mathrm{\πf}}_{d\max }{T}_{s}n\right)]}{1-\frac{1}{N^2}[N+2\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}(N-n)J_0\left({2\mathrm{\πf}}_{d\max }{T}_{s}n\right)]+\frac{{E\left|W\left(m\right)\right|}^2}{{E\left|S\left(m\right)\right|}^2}}$

其中E|S(m)|²、E|W(m)|²分别表示信号及噪声的功率，最大多普勒频偏f_dmax通过f_dmax＝αf_H得到，α为多普勒因子,f_H为信号频率上限。

6.根据权利要求1所述的OFDM水声通信并行迭代ICI消除方法，其特征是：所述的并行迭代MMSE均衡为：

（1）对MMSE准则进行修正

将均衡矩阵H_mmse中的信噪比SNR用信号干扰噪声比SINR代替，即

$H_{\mathrm{mmse}}=R_{H_{\mathrm{ls}}{H}_{\mathrm{ls}}}{(R_{H_{\mathrm{ls}}{H}_{\mathrm{ls}}}+\frac{\mathrm{\β}}{\mathrm{SINR}}I)}^{-1}{H}_{\mathrm{ls}}$

其中H_ls是LS估计得到的信道频响，

是其自相关函数，I是单位矩阵，对于QPSK调制β取1；

（2）进行并行迭代MMSE均衡

ξ为迭代终止阈值，p_max表示最大迭代次数；

①求出在迭代均衡过程中固定不变的矩阵C^r1，C^r2及信号干扰噪声比SINR；

②进行LS估计，获得导频位置处的信道频响

并进行一阶线性频域插值获得所有子载波上的信道频响H_ls，然后借助步骤①中求得的SINR，得到H_mmse；再由

完成MMSE均衡，得到发送信号的估计值

③初始化迭代次数p＝1，令第一次迭代前的接收信号值Y⁽⁰⁾＝Y，Y为接收符号，第一次迭代前的接收信号值估计值

④估计出导频位置处的信道频响

$\hat{H}(m_{\mathrm{ii}},m_{\mathrm{ii}})=\frac{Y^{p-1}\left(m_{\mathrm{ii}}\right)}{P\left(m_{\mathrm{ii}}\right)}=H(m_{\mathrm{ii}},m_{\mathrm{ii}})+\frac{\mathrm{ICI}\left(m_{\mathrm{ii}}\right)+W\left(m_{\mathrm{ii}}\right)}{P\left(m_{\mathrm{ii}}\right)}$

其中P表示导频序列；

⑤获取中间时刻信道冲击响应h^ave，继而得出信道时变斜率α^r1，α^r2；

⑥获取矩阵H_mid，

及

⑦求出矩阵H，H'及H″，并结合前一次估计出的

估计出接收信号中的ICI^(p)分量；

⑧从接收符号Y中将ICI分量减掉，令Y^p＝Y-ICI^(p)；

⑨获得H_mmse再由

完成MMSE均衡，得到p次迭代后的发送信号的估计值

⑩当估计值足够收敛即满足

或p≥p_max时，终止迭代；否则迭代次数p增1，返回到步骤④。

7.根据权利要求6所述的OFDM水声通信并行迭代ICI消除方法，其特征是：p_max小于或等于8。

Images