对于多径传播环境，信道模型通常可以由信号的离开角(AOD，angle of departure)、到达角(AOA，angle of arrival)和传播时延等参数进行表征^[1].已有的研究表明，多径传播环境下的信道估计问题可归结为信道的多维参数估计问题^[2-4]. Bazzi等^[2]在正交频分复用(OFDM，orthogonal frequency division multiplexing)系统中基于二维的多重信号分类算法实现了信号角度和时延的联合估计.然而，多维网格搜索会使该算法产生较高的运算量.针对多输入多输出(MIMO，multiple input multiple output)系统，Adeogun等^[3]提出了一种三维旋转不变子空间(ESPRIT，estimating signal parameters via rotational invariance techniques)算法，实现了信号多普勒频移和角度参数的联合估计. de Almeida等^[4]将MIMO系统的接收信号构造为平行因子^[5-6](PARAFAC，parallel factor)模型，并采用三线性交替最小二乘(TALS，trilinear alternating least square)算法对该模型进行拟合，实现了信号角度和传播时延的联合估计.然而，当所构造的PARAFAC模型的维度较高时，该方案的计算复杂度较高.

利用多维矩阵低秩分解的思想，在多径传播环境下对均匀线阵的接收信号进行建模，构造出具有PARAFAC结构的信号模型，利用截尾奇异值分解(SVD，singular value decomposition)法对该模型进行降维，并采用TALS算法对降维后的信号模型进行拟合，实现了信号角度和传播时延的联合估计.本研究的主要贡献如下：1) 利用截尾SVD法实现了信号模型维度的转换；2) 与传统PARAFAC分解方法相比，所提方法在拟合过程中占用的存储空间更少，拟合速度更快；3) 所提方法无须对角度和传播时延等参数进行配对处理即可分辨出不同路径的参数信息.

1 系统模型

考虑的系统模型如图 1所示，在频率选择性衰落信道环境下，远场用户发送的信号通过P条不同的路径到达配置了阵元数为M_R的均匀线阵.其中，用户设备配置的天线数为M_T=1，接收端阵元间隔为d.系统采用OFDM技术，循环前缀的长度为T_cp，采样间隔为T_S.假设每个OFDM符号块包含N个时隙，其中第一个时隙为导频信号，剩余的N-1个时隙为数据信息.令s(t)=[s₁(t), …, s_K(t)]^T∈$ \mathbb{C}$^K×1为第t个符号块时用户发送的导频向量.其中s_k(t), k=1, 2, …, K，表示第t个符号块时第k子载波上的导频信号.

考虑的信道模型h(t)∈$ \mathbb{C}$^M_R×1为^[4]

$ \mathit{\boldsymbol{h}}\left( t \right) = \sum\limits_{p = 1}^P {{\beta _{t,p}}} \mathit{\boldsymbol{\alpha }}\left( {{\theta _p}} \right)\delta \left( {t - {\tau _p}} \right) $

(1)

其中：β_{t, p}为用户发送第t个符号块时第p条路径的衰减系数；θ_p和τ_p为信号通过第p条路径到达接收阵元时的AOA和传播时延，并且最大传播时延小于T_cp. α(θ)∈$ \mathbb{C}$^M_R×1为接收端的阵列响应矢量，其表达式为

$ \mathit{\boldsymbol{\alpha }}\left( \theta \right) = {\left[ {1, \cdots ,{\text{e}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}\left( {{M_{\rm{R}}} - 1} \right)d\sin \left( \theta \right)/\lambda }}} \right]^{\rm{T}}} $

(2)

其中：λ为信号的波长，d=λ/2.假设AOA和传播时延在连续T₁个OFDM符号块内是保持静止的，且T₁大于用户发送的总符号块数T；衰减系数在每个符号块内是保持恒定的，而在各符号块之间是变化的.因此，用户发送符号块t时，第k个子载波上的信道频域响应表达式为

$ {\mathit{\boldsymbol{h}}_k}\left( t \right) = \sum\limits_{p = 1}^P {{\beta _{t,p}}} \mathit{\boldsymbol{\alpha }}\left( {{\theta _p}} \right){\text{e}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}{\tau _p}{f_{\rm{s}}}\left( {k - 1} \right)/K}} $

(3)

其中f_s为采样频率.去除循环前缀并进行离散傅里叶变换后，接收端在接收符号块t时，所有子载波在导频处的接收信号为

$ \mathit{\boldsymbol{Y}}\left( t \right) = \mathit{\boldsymbol{AH}}\left( t \right){\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{S}}\left( t \right) + {\mathit{\boldsymbol{N}}_1}\left( t \right) $

(4)

其中：A=[α(θ₁), …, α(θ_P)]∈$ \mathbb{C}$^M_R×P是接收端的阵列响应矩阵，B=[b(τ₁), …, b(τ_P)]∈$ \mathbb{C}$^K×P是传播时延的频率响应矩阵，并且有

$ \mathit{\boldsymbol{b}}\left( {{\tau _p}} \right) = {\left[ {1, \cdots ,{\text{e}^{^{ - 2{\rm{ \mathsf{ π} }}{f_{\rm{s}}}{\tau _p}\left( {K - 1} \right)/K}}}} \right]^{\rm{T}}} $

(5)

H(t)=diag(β_{t, 1}, …, β_{t, P})是由第t个符号块内的衰减系数所构成的对角阵；S(t)=diag(s₁(t), …, s_K(t))是由第t符号块的导频向量构成的对角阵；N₁(t)是对应的噪声矩阵.假设所有T个符号块发送的导频向量相同，即S(t)=S，t=1, 2, …, T，令$ \mathit{\boldsymbol{\tilde B}}$=S^TB，则

$ \mathit{\boldsymbol{Y}}\left( t \right) = \mathit{\boldsymbol{AH}}\left( t \right){{\mathit{\boldsymbol{\tilde B}}}^{\rm{T}}} + {\mathit{\boldsymbol{N}}_1}\left( t \right) $

(6)

将T个符号块在导频处的接收信号进行堆叠，则

$ \begin{gathered} \mathit{\boldsymbol{Y}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Y}}^{\rm{T}}}\left( 1 \right)} \\ \vdots \\ {{\mathit{\boldsymbol{Y}}^{\rm{T}}}\left( T \right)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\tilde BH}}\left( 1 \right)} \\ \vdots \\ {\mathit{\boldsymbol{\tilde BH}}\left( T \right)} \end{array}} \right]{\mathit{\boldsymbol{A}}^{\rm{T}}} + {\mathit{\boldsymbol{N}}_2} = \hfill \\ \;\;\;\;\;\;\;\;\;\left( {\mathit{\boldsymbol{H}} \odot \mathit{\boldsymbol{\tilde B}}} \right){\mathit{\boldsymbol{A}}^{\rm{T}}} + {\mathit{\boldsymbol{N}}_2} \in {\mathbb{C}^{TK \times {M_{\rm{R}}}}} \hfill \\ \end{gathered} $

(7)

其中：H∈$ \mathbb{C}$^T×P是T个符号块内所有衰减系数组成的矩阵，有D_t(H)=H(t). D_t(·)表示对角化算子，即选择矩阵的第t行元素作为生成矩阵的主对角线元素，生成矩阵的其余元素为0. ⊙表示矩阵的Khatri-Rao积.式(7) 可视为三维PARAFAC^[5-6]模型$ \mathscr{Y}$∈$ \mathbb{C}$^M_R×T×K沿着时间维度展开的矩阵展开式. $ \mathscr{Y}$的3个加载矩阵分别由A、H和$ \mathit{\boldsymbol{\tilde{B}}}$来表征，即

$ \mathscr{Y} = {\mathscr{T}_p} \times {\;_1}\mathit{\boldsymbol{A}} \times {\;_2}\mathit{\boldsymbol{H}} \times {\;_3}\mathit{\boldsymbol{\tilde B + }}\mathit{\mathscr{N}} $

(8)

其中：$ \mathscr{T}$_P是一个P阶的三维单位阵，$ \mathscr{N}$为对应的三维噪声矩阵.该模型还可写为标量形式

$ y{ _{{m_{\rm{R}}},t,k}} = \sum\limits_{p = 1}^P {{a_{{m_{\rm{R}}},p}}} {\beta _{t,p}}{{\tilde b}_{k,p}} + {n_{{m_{\rm{R}}},t,k}} $

(9)

其中：y_{mR, t, k}和n_{mR, t, k}分别是三维矩阵$ \mathscr{Y}$和$ \mathscr{N}$中的元素，a_{mR, p}、β_{t, p}和$ {{{\mathit{\tilde{b}}}}_{\mathit{k}\rm{, }\mathit{p}}}$分别是加载矩阵A、H和$ {\mathit{\boldsymbol{\tilde{B}}}}$中对应的元素.式(9) 可视为一个三线性分解表达式，表示$ \mathscr{Y}$中的每一个元素均可由P个秩1的三重因子乘积之和表征.在省略噪声项的情况下，$ \mathscr{Y}$还可以沿着空间(接收天线)维度和频率(子载波)维度展开，并有如下2种矩阵展开式,

$ \mathit{\boldsymbol{X}} = \left( {\mathit{\boldsymbol{A}} \odot \mathit{\boldsymbol{H}}} \right){{\mathit{\boldsymbol{\tilde B}}}^{\rm{T}}} \in {\mathbb{C}^{{M_{\rm{R}}}T \times K}} $

(10)

$ \mathit{\boldsymbol{Z}} = \left( {\mathit{\boldsymbol{\tilde B}} \odot \mathit{\boldsymbol{A}}} \right){\mathit{\boldsymbol{H}}^{\rm{T}}} \in {\mathbb{C}^{{M_{\rm{R}}}K \times T}} $

(11)

2 降维PARAFAC分解

TALS算法作为一种无约束的优化算法，由于其能保证单调收敛，简单易行，已被广泛应用于多维矩阵模型的拟合^[5-6].当满足分解唯一性条件时，可采用TALS算法对$ \mathscr{Y}$进行拟合，实现角度和传播时延的联合估计.然而当K≥M_RT＞P时，采用TALS算法对$ \mathscr{Y}$进行拟合的计算复杂度较高.为降低TALS算法的计算复杂度，利用截尾SVD的方法对$ \mathscr{Y}$进行降维，具体实现过程如下：

对X^T进行截尾SVD^[7]，则

$ {\mathit{\boldsymbol{X}}^{\rm{T}}} \approx {\mathit{\boldsymbol{U}}_p}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_p}\mathit{\boldsymbol{V}}_P^{\rm{H}} = {\mathit{\boldsymbol{U}}_p}\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} $

(12)

其中：Σ_P∈$ \mathbb{C}$^P×P是X^T的前P个奇异值构成的对角阵，U_P∈$ \mathbb{C}$^K×P和V_P∈$ \mathbb{C}$^M_RT×P是Σ_P对应的左右奇异值向量构成的酉矩阵.由式(10) 可知，X^T的表达式为

$ {\mathit{\boldsymbol{X}}^{\rm{T}}} = \mathit{\boldsymbol{\tilde B}}{\left( {\mathit{\boldsymbol{A}} \odot \mathit{\boldsymbol{H}}} \right)^{\rm{T}}} $

(13)

当K＞P，并且$ {\mathit{\boldsymbol{\tilde{B}}}}$列满秩时，对比式(12) 与(13) 可知，$ {\mathit{\boldsymbol{\tilde{B}}}}$和U_P的列向量张成的子空间属于同一个子空间，即它们之间存在线性变换的关系.因此，存在一个非奇异转换矩阵Φ∈$ \mathbb{C}$^P×P，使得

$ \mathit{\boldsymbol{\tilde B}} = {\mathit{\boldsymbol{U}}_P}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} $

(14)

令式(13) 两边同时左乘U_P^H，则

$ \mathit{\boldsymbol{U}}_P^{\rm{H}}{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{ = \boldsymbol{\varPhi} }}{\left( {\mathit{\boldsymbol{A}} \odot \mathit{\boldsymbol{H}}} \right)^{\rm{T}}} \in {\mathbb{C}^{P \times {M_{\rm{R}}}T}} $

(15)

此时，$ \mathscr{Y}$转化为一个新的平行因子模$ \mathscr{X}$∈$ \mathbb{C}$^M_R×T×P，$ \mathscr{X}$的加载矩阵分别为A、H和Φ.当K远大于P时，$ \mathscr{X}$中包含的数据量是远小于$ \mathscr{Y}$的. $ \mathscr{X}$可表示为如下3个矩阵展开式：

$ {\mathit{\boldsymbol{X}}_1} = \left( {\mathit{\boldsymbol{H}} \odot \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}} \right){\mathit{\boldsymbol{A}}^{\rm{T}}} $

(16)

$ {\mathit{\boldsymbol{X}}_2} = \left( {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} \odot \mathit{\boldsymbol{A}}} \right){\mathit{\boldsymbol{H}}^{\rm{T}}} $

(17)

$ {\mathit{\boldsymbol{X}}_3} = \left( {\mathit{\boldsymbol{A}} \odot \mathit{\boldsymbol{H}}} \right){\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}} $

(18)

根据式(16)~(18)，可交替最小化下列条件最小二乘(LS，least square)准则来联合估计H、A和Φ：

$ {{\mathit{\boldsymbol{\hat H}}}^{\left( i \right)}} = \mathop {{\rm{arg}}\;{\rm{min}}}\limits_\mathit{\boldsymbol{H}} \left\| {{\mathit{\boldsymbol{X}}_2} - \left( {{{\mathit{\boldsymbol{ \boldsymbol{\hat \varPhi} }}}^{\left( {i - 1} \right)}} \odot {{\mathit{\boldsymbol{\hat A}}}^{\left( {i - 1} \right)}}{\mathit{\boldsymbol{H}}^{\rm{T}}}} \right)} \right\|_F^2 $

(19)

$ {{\mathit{\boldsymbol{ \boldsymbol{\hat \varPhi} }}}^{\left( i \right)}} = \mathop {{\rm{arg}}\;{\rm{min}}}\limits_\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} \left\| {{\mathit{\boldsymbol{X}}_3} - \left( {{{\mathit{\boldsymbol{\hat A}}}^{\left( {i - 1} \right)}} \odot {{\mathit{\boldsymbol{\hat H}}}^{\left( i \right)}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}} \right)} \right\|_F^2 $

(20)

$ {{\mathit{\boldsymbol{\hat A}}}^{\left( i \right)}} = \mathop {{\rm{arg}}\;{\rm{min}}}\limits_\mathit{\boldsymbol{A}} \left\| {{\mathit{\boldsymbol{X}}_1} - \left( {{{\mathit{\boldsymbol{\hat H}}}^{\left( i \right)}} \odot {{\mathit{\boldsymbol{ \boldsymbol{\hat \varPhi} }}}^{\left( i \right)}}{\mathit{\boldsymbol{A}}^{\rm{T}}}} \right)} \right\|_F^2 $

(21)

其中：i表示迭代次数，‖·‖_F表示Frobenius范数，$ {{\mathit{\boldsymbol{\hat H}}}^{\left( \mathit{i} \right)}}{\rm{、}}{{\mathit{\boldsymbol{\hat A}}}^{\left( \mathit{i} \right)}}{\rm{和}}{{\mathit{\boldsymbol{ \boldsymbol{\hat \varPhi} }}}^{\left( \mathit{i} \right)}}$分别表示H、A和Φ的第i次迭代的估计值.根据式(19)~(21)，依次交替拟合上述3个矩阵，即可实现$ \mathscr{X}$的分解.具体过程如下：

步骤1  根据式(12)~(15) 得到$ \mathscr{X}$，初始化矩阵$ {{\mathit{\boldsymbol{\hat A}}}^{\left( {\rm{0}} \right)}}{\rm{、}}{{\mathit{\boldsymbol{ \boldsymbol{\hat \varPhi} }}}^{\left( {\rm{0}} \right)}}$，令i=0；

步骤2  令i=i+1；

步骤3  根据式(19)，计算H的LS估计

$ {{\mathit{\boldsymbol{\hat H}}}^{\left( i \right)}} = {\left( {{{\left( {{{\mathit{\boldsymbol{ \boldsymbol{\hat \varPhi} }}}^{\left( {i - 1} \right)}} \odot {{\mathit{\boldsymbol{\hat A}}}^{\left( {i - 1} \right)}}} \right)}^\dagger }{\mathit{\boldsymbol{X}}}{_2}} \right)^{\rm{T}}}; $

步骤4  根据式(20)，计算Φ的LS估计

$ {{\mathit{\boldsymbol{ \boldsymbol{\hat \varPhi} }}}^{\left( i \right)}} = {\left( {{{\left( {{{\mathit{\boldsymbol{\hat A}}}^{\left( {i - 1} \right)}} \odot {{\mathit{\boldsymbol{\hat H}}}^{\left( i \right)}}} \right)}^\dagger }{\mathit{\boldsymbol{X}}}_{3}} \right)^{\rm{T}}}; $

步骤5  根据式(21)，计算A的LS估计

$ {{\mathit{\boldsymbol{\hat A}}}^{\left( i \right)}} = {\left( {{{\left( {{{\mathit{\boldsymbol{\hat H}}}^{\left( i \right)}} \odot {{\mathit{\boldsymbol{ \boldsymbol{\hat \varPhi} }}}^{\left( i \right)}}} \right)}^\dagger }{\mathit{\boldsymbol{X}}}_{1}} \right)^{\rm{T}}}; $

步骤6  令ϕ⁽ⁱ⁾= $ {\left\| {{\mathit{\boldsymbol{X}}_{\rm{2}}} - \left( {{{\mathit{\boldsymbol{ \boldsymbol{\hat \varPhi} }}}^{\left( \mathit{i} \right)}} \odot {{\mathit{\boldsymbol{\hat A}}}^{\left( \mathit{i} \right)}}} \right){{\mathit{\boldsymbol{\hat H}}}^{\left( \mathit{i} \right){\rm{T}}}}} \right\|_\mathit{F}}$，若|ϕ⁽ⁱ⁾-ϕ^(i-1)|＜ε(ε＜10^-8)，则迭代结束；否则，跳转至步骤2;

步骤7   $ \mathit{\boldsymbol{\hat B}}{\rm{ = }}{\mathit{\boldsymbol{S}}^\dagger }{\mathit{\boldsymbol{U}}_\mathit{P}}\mathit{\boldsymbol{ \boldsymbol{\hat \varPhi} }}$.

其中：ϕ⁽ⁱ⁾表示第i次迭代过程中的代价函数，$ {\mathit{\boldsymbol{ \boldsymbol{\hat \varPhi} }}}$和$ {\mathit{\boldsymbol{\hat B}}}$分别表示Φ和B的估计值，(·)^†为求伪逆运算.对于第p条路径的AOA信息$ {{\mathit{\hat \theta }}_\mathit{p}}$和传播时延$ {{\mathit{\hat \tau }}_\mathit{p}}$，可采用对估计矩阵$ {\mathit{\boldsymbol{\hat A}}}$和$ {\mathit{\boldsymbol{\hat B}}}$的第p列求平均的方法获得.

3 可辨识性分析

令k_A、k_H和k_Φ分别表示加载矩阵A、H和Φ的k秩^[8]，根据PARAFAC模型分解唯一性^[5]可知，若满足

$ {k_\mathit{\boldsymbol{A}}} + {k_\mathit{\boldsymbol{H}}} + {k_\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}} > 2P + 2 $

(22)

则加载矩阵本质唯一.此时，估计结果与加载矩阵之间有如下关系

$ \mathit{\boldsymbol{\hat A}} = \mathit{\boldsymbol{A \boldsymbol{\varPi} }}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_1},\mathit{\boldsymbol{\hat H}} = \mathit{\boldsymbol{H \boldsymbol{\varPi} }}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_2},\mathit{\boldsymbol{ \boldsymbol{\hat \varPhi} }} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} \boldsymbol{\varPi} }}{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}_3} $

(23)

其中：Π是列模糊矩阵；Δ₁、Δ₂和Δ₃是尺度模糊矩阵，并满足Δ₁Δ₂Δ₃=I.对于列模糊和尺度模糊，可采用归一化的方法消除^[5].由式(23) 可知，加载矩阵的估计结果具有相同的列模糊.因此，不同路径的AOA和传播时延能够自动配对.由于A是范德蒙德矩阵，而H为随机矩阵，则A和H具有满k秩.因为S和B是列满秩矩阵，则rank($ {\mathit{\boldsymbol{\tilde{B}}}}$)=rank(S^TB)=rank(B)=P.同理，rank(Φ)=P.又因为矩阵B、S和U_P不存在线性相关列，则Φ具有满k秩.因此，式(22) 可等效为

$ \min \left( {{M_{\rm{R}}},P} \right) + \min \left( {T,P} \right) > P + 2 $

(24)

当系统参数满足式(24) 时，A、Φ和H在LS意义上是可辨识的.

4 仿真结果

使用蒙特卡洛仿真，对角度和时延估计的性能进行分析及验证.假设H为随机生成矩阵，噪声为独立同分布的高斯白噪声.若无特殊说明，假设多径数P=5，到达角θ=(-10°, 20°, 30°, 40°, 60°)，传播时延τ=(0.1T_S, 1T_S, 1.5T_S, 2T_S, 3T_S).导频信号采用正交相移键控调制.信道估计性能由均方根误差(RMSE, root mean squared error)表征.

图 2给出了所提方法对AOA和传播时延的估计效果.其中，M_R=8，T=8，K=64，蒙特卡洛仿真次数Q=30，信噪比为10 dB.由图 2可知，即使在低信噪比情况下，所提方法仍然能够将所有的AOA和传播时延联合估计出来，并且估计结果均集中于原始参数值的附近，具有较好的性能.

图 3给出了所提方法和已有方法性能的比较.已有方法选择了传统PARAFAC分解方法和ESPRIT方法.其中，T=8，K=64，Q=500；当P=3时，θ=(-10°, 20°, 30°)，τ=(0.1T_S, 1T_S, 1.5T_S). 图 3表明，所提方法的估计性能是优于ESPRIT算法的.特别是当接收天线数目M_R小于传播路径个数P时，ESPRIT算法是无法工作的，而所提方法依然具有较好的估计性能，并接近于传统PARAFAC分解方法的估计性能.此外，由于所提方法在进行截断处理时造成了部分信号成分的损失，在低信噪比情况下的估计性能略低于传统PARAFAC分解方法.在中高信噪比情况下两者的RMSE曲线是近乎重合的.

由第2节可知，算法的总复杂度由算法的收敛次数及单次迭代的复杂度所决定.对于传统PARAFAC分解方法，其单次迭代的复杂度为O(P³+M_RTKP)^[5]；而对于所提方法，单次迭代的复杂度为O(P³+M_RTP²).当K远大于P时，后者的运算量是远小于前者的.本研究采用CPU占用时间(单位：s)来表征算法复杂度的具体情况. 表 1给出了2种方法在不同信噪比情况下达到收敛时的CPU占用时间.仿真参数为M_R=8，T=8，K=64. 表 1表明，所提方法具有更快的收敛速度.此外，所提方法在拟合过程中占用的存储空间也更少.

5 结束语

利用多维矩阵低秩分解思想，在采用OFDM技术的阵列天线系统中提出了一种信道多维参数的联合估计方法.借助导频信号，在接收端对接收信号进行建模，构造出三维PARAFAC信号模型，并采用TALS算法对降维后的信号模型进行拟合，实现了信道多个参数的联合估计.计算机仿真及分析验证了所提方法的有效性.