1.   引言

目前，移动通信已进入5G时代，然而当面向无人机、超级高铁、卫星和车联网等超高速移动应用场景(时速超过500 km/h)时^[1]，多普勒频移会对5G标准中所采用的正交频分复用(Orthogonal Frequency Division Multiplexing, OFDM)造成严重的子载波间干扰。于是，研究人员开始关注超高速移动场景下的新型波形^[2]。近年来，正交时频空(Orthogonal Time Frequency Space, OTFS)^[3]调制由于能有效对抗信道的多普勒效应而逐渐成为学术界的研究热点^[4]，它继承了OFDM^[5]系统的抗多径能力，并且作为一种基于时延-多普勒(Delay-Doppler, DD)域的新型调制方式，还能使得每个符号即使在高多普勒频移条件下也能经历同样的信道增益，因此大大提高了分集增益^[6]。当数据被承载在资源网格的DD域中^[7]，接收信号可以表示为DD域信道与OTFS传输符号的2维卷积形式。然而，当同时存在多普勒频移和多径效应时，2维卷积形式会导致符号在资源网格上沿时延和多普勒频移两个方向上进行扩散，即引起了多普勒间干扰(Inter-Doppler Interference, IDI)^[8]，这给OTFS接收机的设计带来极大的挑战。

为了克服IDI问题，在信道估计方面，文献[9,10]都采用了一种单脉冲方法，并在该单脉冲周围添加保护间隔来区分导频和数据引起的干扰，然而该算法估计的是信道采样后分布在DD域的脉冲形式，该方案并不能解决IDI带来能量泄露问题。文献[11]设计了一组变换域基函数，通过构造一个低维子空间来建模OTFS信道，能很好地估计信道每符号内采样点的变化，但是文中采用线性最小均方误差(Linear Minimum Mean Square Error, LMMSE)算法来进行低维子空间的信道估计，这需要信道的先验信息，在实际场景中很难获取。信道均衡也是OTFS接收机设计的重要功能，也是提升解调性能的关键。文献[12]利用OTFS信道矩阵的块循环特性参与解调过程，从而提出了低复杂度线性接收机。此外，Thaj等人^[13]分析了OTFS传输矩阵在不同域的结构，并基于此提出了基于最大比合并(Maximal Ratio Combining, MRC)的迭代检测方案，该方案能有效获取信道的分集增益。这两种方案都能在一定程度上取得良好的均衡效果，然而这两种方案都适用于单一散射体环境，对于实际中存在的多散射体环境，即多普勒频移不断变化的情况并不适用。

针对以上问题，本文从提升接收机的整体性能出发，设计了一种联合信道估计和信道均衡的接收机，并提出一种基系数到信道等效矩阵的低复杂度转换算法，在降低接收机复杂度的同时能取得良好的通信性能。本文的主要贡献如下：

(1)针对OTFS接收机面临的IDI问题，将OTFS信道建模成能反映多普勒域采样点变化的DD域信道响应，然后用基扩展模型(Basis Expansion Model, BEM)将信道估计问题转换成基系数恢复问题。

(2)进一步提出稀疏贝叶斯(Sparse Bayesian Learning, SBL)方法来估计信道的基系数和信道噪声，以更准确地估计包括多普勒采样点在内的DD域信道和更好地消除IDI的干扰。

(3)为了对DD域接收符号进行均衡，提出一种高效的基系数转换方法，并进行了理论推导，通过卷积和的形式实现基系数到DD域等效矩阵的获取。

(4)针对SBL估计得到的信道噪声，设计了一种基于阻尼最小二乘最小残差(damped Least Square Minimum Residual, d-LSMR)的信道均衡器，该均衡器能在保证低复杂度计算的前提下进一步消除信道噪声以及IDI干扰。

(5)采用3GPP指定的扩展车辆信道模型(Extended Vehicular A model, EVA)，充分模拟了富含多散射体且多普勒频移不断变化的高速移动场景，并通过系统仿真结果验证了所提方案的优越性。

2.   OTFS系统模型

2.1   传输模型

OTFS将数据流映射到DD域进行数据传输，图1展示了OTFS系统的传输模型。在发射端，数据流经过调制之后与导频序列一起被映射到DD域的资源网格${\boldsymbol{X}} \in {\mathbb{C}^{M \times N}}$上，这里M,N分别表示时延域和多普勒域的资源元素数量。经过辛有限傅里叶逆变换(Inverse Symplectic Finite Fourier Transform, ISFFT)之后，即${\boldsymbol{D}} = {{\boldsymbol{F}}_M}{\boldsymbol{XF}}_N^{\text{H}}$，${\boldsymbol{X}}$从DD域转换到了TF域，${\boldsymbol{D}} \in {\mathbb{C}^{M \times N}}$代表了TF域的资源网格，${{\boldsymbol{F}}_M} \in {\mathbb{C}^{M \times M}}$和${{\boldsymbol{F}}_N} \in {\mathbb{C}^{N \times N}}$是归一化快速傅里叶变换矩阵。接着，与OFDM调制操作类似，本文需要对${\boldsymbol{D}}$进行IFFT操作将其变换到时域，即${\boldsymbol{S}} = {\boldsymbol{F}}_M^{\text{H}}{\boldsymbol{D}} = {\boldsymbol{XF}}_N^{\text{H}}$，这里${\boldsymbol{S}} \in {\mathbb{C}^{M \times N}}$。为了让发送的数据能在信道中进行传输，需要将${\boldsymbol{S}}$转换成串流数据，即${\boldsymbol{s}} = {\text{vec(}}{\boldsymbol{S}}{\text{) = (}}{\boldsymbol{F}}_N^{\text{H}} \otimes {{\boldsymbol{I}}_M}{\text{)}}{\boldsymbol{x}}$，这里${\boldsymbol{s}} \in {\mathbb{C}^{MN \times 1}}$, ${\boldsymbol{x}} = {\text{vec}}\left( {{{\boldsymbol{X}}^{\text{T}}}} \right) \in {\mathbb{C}^{MN \times 1}}$代表了DD域的串流发射信号，${{\boldsymbol{I}}_M} \in {\mathbb{C}^{M \times M}}$代表了单位矩阵。

图 1  OTFS传输系统

下载: 全尺寸图片 幻灯片

综上，可以得到OTFS系统中DD域的输入输出关系

其中，${\boldsymbol{H}} = \left( {{{\boldsymbol{I}}_M} \otimes {{\boldsymbol{F}}_N}} \right){{\boldsymbol{P}}^{\text{H}}}{\boldsymbol{GP}}\left( {{{\boldsymbol{I}}_M} \otimes {\boldsymbol{F}}_N^{\text{H}}} \right) \in {\mathbb{C}^{MN \times MN}}$是DD域的等效信道矩阵，${\boldsymbol{P}} \in {\mathbb{C}^{MN \times MN}}$是置换矩阵^[13]，${\boldsymbol{\tilde w}} = \left( {{{\boldsymbol{F}}_N} \otimes {{\boldsymbol{I}}_M}} \right){\boldsymbol{w}} \in {\mathbb{C}^{MN \times 1}}$代表DD域的信道噪声。

2.2   信道的基扩展模型

在时域中，输入信号和时域信道进行卷积运算得到输出信号；类似地，在DD域中，输出信号是由输入信号与DD域信道的2维卷积得到，将其表示为矩阵形式为

其中，${\boldsymbol{y}} = {\text{vec}}({{\boldsymbol{Y}}^{\text{T}}}) \in {\mathbb{C}^{MN \times 1}}$表示DD域的串行接收信号，${\boldsymbol{h}} \in {\mathbb{C}^{MNL \times 1}}$代表DD域信道响应的向量形式，$\mathcal{X} = [{\mathcal{X}_0}, \cdots ,{\mathcal{X}_l}, \cdots ,{\mathcal{X}_{L - 1}}] \in {\mathbb{C}^{MN \times MNL}}$是由DD域的发射向量组成的，其中的$\mathcal{X}_1 \in {\mathbb{C}^{MN×MN}}$是一个块对角矩阵，包含了$M$个维度为$N \times N$的小矩阵${{\boldsymbol{\varDelta }}_i}$，其中 ${{\boldsymbol{\varDelta }}_i}$是由发射向量组成的循环矩阵。

从式(2)可以看出，待求解参数${\boldsymbol{h}}$的维度大于观测值${\boldsymbol{y}}$的维度，这说明该方程是一个欠定方程，不能直接求解。为了估计${\boldsymbol{h}}$，本文利用BEM模型来降低待求解的参数量的个数。BEM可以利用一组固定的基函数和其基系数的线性组合来拟合信道的时域响应，由此将信道系数的估计转变为基系数的估计。在本文中，离散长椭球序列基扩展模型(Discrete Prolate Spheroidal Basis Expansion Model, DPS-BEM)^[14]被用来对时域信道抽头进行建模

其中，${{\boldsymbol{b}}_q} \in {\mathbb{C}^{MN \times 1}}$为基向量，${{\boldsymbol{c}}_l} \in {\mathbb{C}^{Q\times 1}}$表示基系数，$Q$为基向量的维度，${{\boldsymbol{\varepsilon }}_l} \in {\mathbb{C}^{MN\times 1}}$表示误差向量。根据时域信道和DD域信道的关系，将时域信道抽头，转换到DD域可以得到如式(4)的矩阵表达式

将式(4)代入式(2)中可以得到观测信号与基系数的关系为

其中，${\boldsymbol{c}} = {[{{\boldsymbol{c}}_0}, \cdots ,{{\boldsymbol{c}}_l}, \cdots ,{{\boldsymbol{c}}_{L - 1}}]^{\text{T}}} \in {\mathbb{C}^{QL\times 1}}$代表基系数矩阵，${\boldsymbol{\varPhi }} \triangleq \mathcal{X}({{\boldsymbol{I}}_L} \otimes {\boldsymbol{A}}) \in {\mathbb{C}^{MN \times QL}}$。

3.   基于稀疏贝叶斯的基系数和噪声估计

2.2节中，BEM将信道估计问题转换成基系数估计问题。考虑到在实际进行信道估计时路径未知，待估计的路径数远远大于真实路径，因此待估计的基系数也呈稀疏性。为了利用基系数的稀疏性，同时能准确地估计噪声，本节提出基于SBL的基系数估计方法来对信道和噪声同时进行估计。

在式(5)中，已知${\boldsymbol{\tilde w}}$服从$\mathcal{N}({0},{\sigma ^2}{\boldsymbol{I}})$分布，则 $p\left( {{\boldsymbol{y}}\left| {{\boldsymbol{c,\tilde w}}} \right.} \right) = \mathcal{N}\left( {{\boldsymbol{\varPhi c}},{\sigma ^2}{\boldsymbol{I}}} \right)$。由于模型中的求解参数和训练参数相当，直接求解可能会造成过拟合，因此我们对待估计参数${\boldsymbol{c}}$做出先验假设

其中，${\boldsymbol{\varLambda}} = {\text{diag(}}[\alpha (0),\alpha \left( 1 \right), \cdots ,\alpha (QL - 1)]{\text{)}} \in {\mathbb{C}^{QL \times QL}}$, ${\boldsymbol{\alpha }} \in {\mathbb{C}^{QL}}$服从伽马分布，则可以利用最大后验概率准则来求得${\boldsymbol{c}}$的最大后验概率服从高斯分布^[15]，其均值和方差分别为

对于未知的超参数${\boldsymbol{\alpha }}$, ${\sigma ^2}$，可以利用Type-Ⅱ 最大似然算法求解，并得到两个超参数的更新准则为

其中，${N_{\rm{p}}}$代表待估计的导频个数，${\varSigma _{jj}}$代表${{\boldsymbol{\varSigma }}_{\boldsymbol{c}}}$对角线元素上的第$j$个值。$\bar {\boldsymbol{\alpha}}$在经过${I_s}$次迭代后收敛。为了便于表示，记$\bar {\boldsymbol{\alpha}} = {{\boldsymbol{\alpha}} ^{ - 1}}$, $\beta = {\sigma ^2}$，基于SBL的基系数估计算法流程见算法1。

算法1　基于SBL学习的基系数估计
输入：${\boldsymbol{y}}$,${\boldsymbol{\varPhi } }$
输出：${\boldsymbol{\tilde c}}$, $\beta$
(1) 初始化：${\bar {\boldsymbol{\alpha} } ^{(0)} } = {{\text{0}}}$, $\beta = 1$
(2) while $||{\bar {\boldsymbol{\alpha} } ^{(i)} } - {\bar {\boldsymbol{\alpha} } ^{(i - 1)} }|| < {\varepsilon _{\boldsymbol{\alpha}} }$
(3) 　for $j = 1,2, \cdots ,{{QL} }$
(4) 　　超参数更新：${\bar \alpha ^{\left( i \right)} }\left( j \right) = \varSigma _{jj}^{(i - 1)} + {\left( { {\boldsymbol{\mu} } _{\boldsymbol{c}}^{\left( {i - 1} \right)}\left( j \right)} \right)^2}$
(5) 　end
(6) 　噪声更新：${\beta ^{(i)} } = \frac{ {||{\boldsymbol{y} } - {\boldsymbol{\varPhi \mu } }_c^{(i)}||_2^2} }{ { {N_{\boldsymbol{p}}} - {\boldsymbol{\varSigma } }_{\boldsymbol{c}}^{(i)}(j)(1 - {\alpha ^{(i)} }(j)\varSigma _{jj}^{(i)})} }$
(7) 　后验更新：通过式(7)和式(8)，更新后验参数${\boldsymbol{\varSigma}} _{\boldsymbol{c}}^{\left( i \right)}$和${\boldsymbol{\mu}} _{\boldsymbol{c}}^{\left( i \right)}$
(8) 　$i = i + 1$
(9) end

下载: 导出CSV 

| 显示表格

4.   基系数转换

在时域上，发送符号表示为接收符号与信道的卷积关系，而基扩展模型对每条信道抽头进行压缩，表示为式(3)的形式。此时，接收信号${\boldsymbol{r}}$可以表示为

其中，${\bar {\boldsymbol{\varXi}} _q} = {\text{diag}}\left( {{{\boldsymbol{b}}_q}} \right) \in {\mathbb{C}^{N \times N}}$为对角矩阵，${{\boldsymbol{C}}_q} \in {\mathbb{C}^{N \times N}}$类似于${\boldsymbol{G}}$，是由${{\boldsymbol{c}}_q}$组成的循环。${\boldsymbol{G}}$矩阵可以通过$\displaystyle\sum\nolimits_{q = 0}^{Q - 1} {{{\bar {\boldsymbol{\varXi}} }_q}{{\boldsymbol{C}}_q}}$得到，这种卷积-乘积和的形式能够通过快速向量运算得到${\boldsymbol{G}}$矩阵。根据DD域信号和时域信号的关系，可以得到DD域接收符号${\boldsymbol{y}}$为

对式(12)中和${\boldsymbol{H}}$有关的项进行分析

其中，${{\boldsymbol{\varXi}} _q}$为块对角矩阵，其中对角矩阵都是一个独立的循环矩阵，${\bar {\boldsymbol{C}}_q}$为块循环矩阵，且每个块矩阵都为一个对角矩阵，将${\bar {\boldsymbol{C}}_q}$表示为

${{\boldsymbol{c}}_q} \in {\mathbb{C}^L}$为基系数。最后代入式(13)中

其中，$\odot$表示Hadamard积。

分析式(15)${\boldsymbol{H}}$的组成，包括两个重要的部分：向量点乘和矩阵加法。上述步骤均没有涉及矩阵乘法运算，因此该方法较为节省空间、高效，能够快速得到频域矩阵，为后续频域均衡算法的优化奠定了基础。

5.   低复杂度迭代均衡器

在数学上，一些著名的迭代方法，如LSMR^[16]、最小残差^[17]和最小二乘QR分解(Least-Squares QR factorization, LSQR)^[18]，被发展来解决大规模的和稀疏的LS问题。该方法采用递归策略，在krylov子空间迭代获得LS解，其精度与LMMSE方法相似。同时，由于${\boldsymbol{H}}$的稀疏性，大大降低了迭代过程的计算量。LSQR算法通过求解式(16)的鞍点问题来求解线性方程组

其中，${\boldsymbol{ \mathcal{H}}}$是对称矩阵。根据Lanczos关系，展开$\mathcal{H}$，可以得到${\boldsymbol{H}}$, ${{\boldsymbol{\nu }}_i}$和${{\boldsymbol{\omega }}_i}$的递推式

其中，${{\eta }_i}$和${{\xi }_i}$为归一化因子，分别使得$||{{\boldsymbol{\nu }}_i}|{|_2} = 1$和$||{{\boldsymbol{\omega }}_i}|{|_2} = 1$，且${{\eta }_1}{{\boldsymbol{\nu }}_1} = {\boldsymbol{y}}$, ${{\xi }_1}{{\boldsymbol{\omega }}_1} = {{\boldsymbol{H}}^{\text{H}}}{{\boldsymbol{\nu }}_1}$。记${\boldsymbol{\mathcal{V}}}{\text{}}_{m+1}= \left[{{\boldsymbol{\nu}} }_{1},{{\boldsymbol{\nu }}}_{2},\cdots ,{{\boldsymbol{\nu}} }_{m+1}\right]$是Krylov子空间中的一组基，${\boldsymbol{\mathcal{W}}}{\text{}}_{m}=\left[{{\boldsymbol{\omega}} }_{1},{{\boldsymbol{\omega}} }_{2},\cdots ,{{\boldsymbol{\omega}} }_{m}\right]$是约束空间中的一组基，${\boldsymbol{H}}$与${\boldsymbol{\mathcal{V}}}{\text{}}_{m+1}$, ${\boldsymbol{\mathcal{W}}}{m}_{}$的关系为

其中，${{\boldsymbol{B}}_{m + 1}} \in {\mathbb{C}^{\left( {m + 1} \right) \times m}}$为双对角阵，式(19)和式(20)是Golub-Kahan双对角化过程^[19]。类似地，LSMR也采用Golub-Kahan双对角化过程产生两组基向量，同时为了考虑消除噪声干扰，可以添加阻尼参数来达到正则化目的，d-LSMR能够考虑到噪声的影响，此时问题转换为^[16]

其中，${{\boldsymbol{\bar R}}_{m + 1}}$通过QR分解得到

算法2描述了利用d-LSMR进行迭代反馈均衡的过程。

算法2　基于d-LSMR的基系数均衡
输入：${\boldsymbol{y}}$,${\boldsymbol{\tilde c}}$, ${\sigma ^2}$
输出：${\boldsymbol{\tilde x}}$
(1) for $q = 1,2, \cdots ,Q$
(2) 　　${\boldsymbol{H = H} }{ + }\left( { {\text{circ} }\left( {\left( { {{\boldsymbol{\varXi}} _q} \odot {c_q}\left( 0 \right)} \right), \cdots \left( { {{\boldsymbol{\varXi}} _q} \odot {c_q}\left( {L - 1} \right)} \right)} \right)} \right)$
(3) end
(4) 初始化：${{\eta }_1}{{\boldsymbol{\nu }}_1} = {\boldsymbol{y}}$, ${{\xi }_1}{{\boldsymbol{\omega }}_1} = {{\boldsymbol{H}}^{\text{H}}}{{\boldsymbol{\nu }}_1}$, $i = 0$ (d-LSMR均衡)
(5) while $i < M$ or ${\varepsilon _i}$收敛
(6) 　　${{\eta }_{i + 1}}{{\boldsymbol{\nu }}_{i + 1}} = {\boldsymbol{H}}{{\boldsymbol{\omega }}_i} - {{\xi }_i}{{\boldsymbol{\nu }}_i}$
(7) 　　${{\xi }_{i + 1}}{{\boldsymbol{\omega }}_{i + 1}} = {{\boldsymbol{H}}^{\text{H}}}{{\boldsymbol{\nu }}_{i + 1}} - {{\eta }_{i + 1}}{{\boldsymbol{\omega }}_i}$
(8) 　　根据式(22)QR分解得到${{\boldsymbol{\bar R}}_{m + 1}}$
(9) 　　求解(21)得到${{\boldsymbol{z}}_i}$，并计算残量
(10) 　 $i = i + 1$
(11) end
(12) $\tilde{{\boldsymbol{x}}}={\boldsymbol{\mathcal{W}}}{\text{} }_{i-1}{{\boldsymbol{z}}}_{i-1}$

下载: 导出CSV 

| 显示表格

6.   复杂度和仿真分析

6.1   复杂度

本节对经典的信道估计和均衡算法的复杂度进行了分析，对比分析如表1所示。

表 1  信道估计/均衡算法复杂度对比

	算法	时间复杂度
信道估计	BEM-LS	$O\left( {{M^2}{N^2}\left( {QL} \right)} \right)$
	BEM-LMMSE	$O\left( {{M^2}{N^2}\left( {QL} \right)} \right)$
	BEM-SBL	$O\left( {{I_S}{M^2}{N^2}\left( {QL} \right)} \right)$
信道均衡	ZF	$O\left( {{M^3}{N^3}} \right)$
	LMMSE	$O\left( {{M^3}{N^3}} \right)$
	d-LMSR	$O\left( {2{M^2}{N^2} + I_K^3} \right)$

下载: 导出CSV 

| 显示表格

6.2   仿真分析

本节通过系统仿真，以此来评估所提信道估计与均衡算法的性能。表2给出了系统仿真的参数设置。表2给出了系统仿真的参数设置。

表 2  仿真系统参数

参数名称	参数值
子载波个数($M$)	32
符号个数($N$)	16
载波频率	5.9 GHz
调制方式	QPSK/16QAM
用户移动速度	121.5～607.5 km/h
子载波间隔	15 kHz
信道模型	EVA[20]
CP长度	7

下载: 导出CSV 

| 显示表格

图2和图3分别呈现了SBL和d-LSMR算法在信噪比为15 dB，归一化多普勒频率为3%和15%时，残量范数随着迭代过程的变化示意图。从图2和图3可以看出，在高速和低速场景下，这两种算法均能较快收敛，其中SBL在10次以内达到收敛，d-LSMR算法在多普勒频移较高时，收敛较慢，但也能在30次以内达到收敛。对于OTFS阶数$MN$而言，30次仍然远远小于$MN$，因此这两种算法的复杂度就如表中分析的一样，并没有因为迭代次数就增加算法的复杂度。

图 2  SBL算法的收敛过程

下载: 全尺寸图片 幻灯片

图 3  d-LSMR算法的收敛过程

下载: 全尺寸图片 幻灯片

从图4可以看出，在低速时，由于信道受到多普勒频移影响较小，变化缓慢，所有的插值算法均能在该场景下达到良好的性能。在高速时，基于完美信道估计的最大比合并(Perfect-Maximal Ratio Combining, P-MRC)和基于完美信道估计的最小均方误差算法(Perfect-Minimum Mean Square Error, P-MMSE)由于使用的是完美的信道估计，因此相比SBL-dLSMR接收机，整体规避了信道估计带来的误差。然而，在这种情况下，SBL-dLSMR接收机的误码率性能却好于P-MRC，并与P-MMSE接近，这说明所提出的信道估计模型用SBL取得了良好的信道估计性能，并且通过基系数转换之后，所采用的dLSMR均衡算法进一步消除了IDI的影响，使整体接收机的性能接近完美信道估计加MMSE均衡方式。从图5中可以看到，本文所提算法在16QAM调制下依然能保持良好误码率性能同时，在两种调制方式下，随着多普勒频移增大，在EVA这种多散射体场景中，P-MRC接收机的性能出现了下降，而本文所提接收机仍然能保持良好的性能趋势，这说明本文所设计的接收机的有效性和应对高速场景的鲁棒性。

图 4  调制方式为QPSK时接收机BER性能

下载: 全尺寸图片 幻灯片

图 5  调制方式为16QAM时接收机BER性能

下载: 全尺寸图片 幻灯片

7.   结论

针对OTFS信号在双选衰落信道中传输造成的IDI问题，本文设计了一个完整的OTFS接收机，重点研究了信道估计、信道均衡以及其中信道模型的转换关系。首先分析了OTFS的传输模型，在此基础上为了能更准确地进行信道估计，本文将DD域信道响应用BEM来进行降维表征，并进一步用SBL来对信道和噪声进行估计。在基系数估计完成之后，进一步提出了基系数和DD域等效信道的低复杂度转换算法，完成了信道估计和均衡的模型衔接。在此基础上，引入基于d-LSMR的信道均衡器，利用了信道估计中的噪声来进一步消除干扰，实现低复杂度的信道均衡。仿真结果验证了所提出的OTFS接收机在富含散射体的高速移动场景下能够理想地解调16QAM调制信号，证明了所提接收机的有效性和鲁棒性。