CN107645460A - 实值平行因子分解的多径参数估算方法 - Google Patents

Abstract

本发明公开了一种实值平行因子分解的多径参数估算方法，其首先构建阵列数据的PARAFAC模型，然后利用前后平滑技术和酉变换技术对接收数据进行处理，构造阵列的增广输出的PARAFAC模型。再利用交替最小二乘算法获得相关导引矢量的估计，最后再利用阵列的旋转不变特性恢复源信号的角度与时延参数。所提算法在进行交替最小二乘时仅涉及实数运算，故相比已有复数算法，本发明算法其计算复杂度更低，且在低信噪比或者小快拍时能获得更优的估计性能，同时还能自动配对所估计角度和时延参数，并有效应对相干源场景，不需要对接收数据进行奇异值分解或者是谱峰搜索运算，在实际的无线通信场景中具有良好的应用前景。

Description

实值平行因子分解的多径参数估算方法

技术领域

本发明涉及一种阵列信号处理技术，具体的说涉及一种实值平行因子分解的多径参数估算方法。

背景技术

伴随着移动通信领域的不断成功与发展，源定位技术的市场需求愈加强烈，从而引起国内外学者们的广泛兴趣。源定位技术是众多无线通信应用领域的基础技术，如无线定位，雷达探测，智能交通和生物医学等。无线定位技术经常涉及到联合频率、角度和时延等参数的估计，本发明主要讨论联合角度和时延估计(joint angle and delayestimation，JADE)。典型的JADE算法主要有多重信号分类(multiple signalclassification，MUSIC)算法，加权子空间拟合(weighted subspace fitting，WSF)算法、基于旋转不变技术的参数估计(estimation signal parameters via the rotationalinvariance techniques，ESPRIT)算法。然而，MUSIC算法需要进行高维谱峰搜索，因而计非常大；WSF算法需要利用最大似然估计，其迭代的复杂度也非常高。且MUSIC算法和WSF算法均存在离散网格不匹配的问题，其估计精度往往有限。ESPRIT算法能充分利用接收信号的空-时信息，然而其存在过载问题(目标数大于阵元数)。平行因子分析(parallel factoranalysis，PARAFAC)，又被称为平行因子分解或者三线性分解，是一种高效的三维低秩张量分解算法，在阵列信号处理中被广泛应用。PARAFAC算法采用迭代的方式进行参数估计，可被视为ESPRIT算法的更一般的表达形式，相比ESPRIT算法，其计算精度往往较高，计算复杂度较小。特别是在高维数据计算时，这种优势更显著。文献张小飞,徐大专.一种新的盲联合角度-时延估计方法[J].哈尔滨工业大学学报,2006,38(11):1893-1897.中提出将PARAFAC算法应用于JADE，取得了优于ESPRIT的估计效果。为降低大规模阵列的计算复杂度，同时利用阵列多维结构，文献F.Q.Wen,G.Zhang.D.Ben.Estimation of multipath parametersin wireless communications using multi-way compressive sensing[J].Journal ofSystems Engineering and Electronics,2015,26(5):908-915中提出基于三维压缩感知的JADE算法。尽管其计算复杂度较小，其估计精度会明显劣于PARAFAC算法。此外，上述两篇文献提出的算法均涉及复数运算，一般来说，复数乘法的运算量是实数乘法运算量的四倍。

发明内容

鉴于以上原因，有必要提供一种将PARAFAC算法扩展为实数运算，实现大大降低参数估计复杂度的实值平行因子分解的多径参数估算方法。

本发明提供一种实值平行因子分解的多径参数估算方法，所述实值平行因子分解的多径参数估算方法包括如下步骤：

S1、构建接收数据的过采样矩阵，根据过采样矩阵构建出过采样矩阵阵列数据的PARAFAC模型；

S2、利用前后平滑技术和酉变换技术对过采样矩阵阵列数据的PARAFAC模型进行处理，构造过采样矩阵阵列的实数增广输出的PARAFAC模型；

S3、通过交替最小二乘算法获得实数增广输出的PARAFAC模型的导引矢量的估计；

S4、根据阵列的旋转不变特性，通过所述导引矢量的估计恢复源信号的角度与时延参数。

本发明提出一种实值平行因子分解的多径参数估算方法。首先构建阵列数据的PARAFAC模型，然后利用前后平滑技术和酉变换技术对接收数据进行处理，构造阵列的增广输出的PARAFAC模型。再利用交替最小二乘算法获得相关导引矢量的估计，最后再利用阵列的旋转不变特性恢复源信号的角度与时延参数。所提算法在进行交替最小二乘时仅涉及实数运算，故相比已有复数算法，本发明算法的计算复杂度更低，且由于前后平滑技术会扩展阵列接收数据，因而所提算法在低信噪比或者小快拍时能获得更优的估计性能，同时还能自动配对所估计角度和时延参数，并有效应对相干源场景，不需要对接收数据进行奇异值分解或者是谱峰搜索运算，在实际的无线通信场景中具有良好的应用前景。

附图说明

图1是非相干源背景下本发明所述实值平行因子分解的多径参数估算方法估计的散点图；

图2是相干源背景下本发明所述实值平行因子分解的多径参数估算方法估计的散点图；

图3是相干源条件下本发明所述实值平行因子分解的多径参数估算方法与与其他算法的角度估计RMSE性能比较示意图；

图4是非相干源条件下本发明所述实值平行因子分解的多径参数估算方法与与其他算法的时延估计RMSE性能比较示意图；

图5是相干源条件下本发明所述实值平行因子分解的多径参数估算方法与与其他算法的角度估计RMSE性能比较示意图；

图6是相干源条件下本发明所述实值平行因子分解的多径参数估算方法与与其他算法的时延估计RMSE性能比较示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明，应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明提供一种实值平行因子分解的多径参数估算方法，所述实值平行因子分解的多径参数估算方法包括如下步骤：

S1、构建接收数据的过采样矩阵，根据过采样矩阵构建出过采样矩阵阵列数据的PARAFAC模型；

S2、利用前后平滑技术和酉变换技术对过采样矩阵阵列数据的PARAFAC模型进行处理，构造过采样矩阵阵列的实数增广输出的PARAFAC模型；

S3、通过交替最小二乘算法获得实数增广输出的PARAFAC模型的导引矢量的估计；

S4、根据阵列的旋转不变特性，通过所述导引矢量的估计恢复源信号的角度与时延参数。

具体的，一个张量就是一个多维向量，为方便读者阅读及理解本发明算法，首先引入关于张量操作的几个基本的定义：

定义1(张量展开)：令为一个N阶张量，的模-n(n＝1,…,N)矩阵展开表示为其中，位于张量的(i₁,…,i_n)位置的元素成为位于矩阵的(i_n,j)处的元素，且

定义2(模-n张量与矩阵乘积)：定义N阶张量与矩阵的模-n乘积为其中且

定义3(张量模乘性质)：N阶张量的模乘性质主要有如下两条：

式中，(·)^T表示转置。

假设通过一个信道传输数字序列s_l，同时利用M根均匀线性阵列的接收天线测量该信源。不考虑阵列接收噪声，则在t时刻天线的响应为

其中，T为符号周期,h(t)是如下信道响应

上式中K为多径数目，θ_k，β_k及τ_k分为第k条衰落路径所对应的波达角(direction-of-arrival，DOA)，权值系数和时延。g(t)是一个已知的脉冲形成函数。a(θ_k)＝[1,exp{-j2πd sinθ_k/λ},…,exp{-j2π(M-1)dsinθ_k/λ}]^T为接收导引矢量，d为天线间距，λ为发射载波波长。假设信道具有有限的时间支撑t∈[0,LT),则最长信道响应为LT＝L_gT+τ_max，L_gT与τ_max分别为脉冲形成函数的时宽和最大多径时延。假设接收信号x(t)以采样率D倍符号率被过采样，且接收信号完全占据N个符号的时长，则可构造一个MD×N的样本矩阵：

定义方向矩阵A＝[a(θ₁),…,a(θ_K)]∈C^M×K,显然，它是一个范德蒙(Vandermonde)矩阵。系数矩阵B＝diag(β₁,…,β_K)∈C^K×K，diag(·)表示对角化操作。令C∈C^K×N为编码矩阵，信源矩阵S＝(BC)^T∈C^N×K。将已知波形g(t)的样本堆叠成一个矩阵g₀＝[g(0),g(1/D),…,g(L-1/D)]，定义波形时延矩阵为G＝[g(τ₁),…,g(τ_K)]∈C^LD×K，g(τ_k)＝[g(0-τ_k),g(1/D-τ_k),…,g(L-1/D-τ_k)]^T是包含g(t-τ_k)波形的所有样本。根据背景技术中两篇文献所提出的JADE算法，表达式(5)可以被表述成

其中⊙表示Khatri-Rao积。

时延估计的核心是基于傅里叶变换将时移转变为相位。g₀的LD点离散傅里叶变换(discrete Fourier transform，DFT)可以被表示为

其中φ＝exp(-j2π/DL)。由DFT的性质可知因此的DFT输出为

式中F即为时延矩阵

上式中φ_k＝exp(-j2πτ_k/L)(k＝1,2…,K)。可以看出，F也是一个范德蒙矩阵.由于是一个已知的对角矩阵，为表达简便，可将其从中分离。利用Tucker张量模型，可以被表述为一个秩为K的三阶张量模型

上式中，表述维数为K×K×K的单位张量。结合定义1可以看出即式(8)中的矩阵为式(10)中的张量的模-3展开。实际上，式(8)即为PARAFAC模型的矩阵表示形式，式(10)为PARAFAC模型的张量表示形式，二者在一定程度上是等价的。

其中，式(10)中的张量为复数张量，如果直接对其进行分解，计算量较大。考虑到A、F的范德蒙特性，可以利用前后平滑技术和酉变换技术对张量数据进行实数化处理，降低参数估计的复杂度。显然，矩阵A与F存在如下特性

式中，Π_M表示反向交换矩阵(下标表示矩阵的维数)，其反对角线元素为1，其余元素为0；令

结合式(11)和定义3中的性质有

令∪_n表示将两个数据张量按照模-n方向进行堆叠，则通过下式可构造一个扩展的PARAFAC模型

可以看出，扩展的PARAFAC模型具有中心对称结构。此时，通过酉变换技术，可以将上述复数张量变换为实值张量

其中，酉变换矩阵具体形式如下：

以M为基数和导引矢量a(θ_k)为例，酉变换后的导引矢量为

因此通过酉变换可将一个复数向量变换为一个实数向量。结合表达式(14)、表达式(15)和定义3有

上式中，此时，是一个维数为M×LD×2N的张量，相比其维数扩大了一倍。此外，扩展的数据仍然为一个PARAFAC模型。

由式(17)可知，构成的方向矩阵和中分别包含了目标的角度与时延信息。通过对进行张量分解即可获得对这些矩阵的估计，从而可以获得对相关参数的估计。本发明采用交替最小二乘算法对进行分解，具体过程如下节所述。

根据定义1，表达式(17)中的张量可以展开成如下矩阵的形式：

上式即为三线性分解模型的矩阵表达形式。Z₁、Z₂和Z₃分别可被视为将张量数据沿着信源方向、时域方向和空域方向展开而获得的矩阵。

交替最小二乘算法是一种高效的三线性分解模算法，其采用最小二乘(leastsquares，LS)代价函数依次交替的拟合三个矩阵，直到拟合误差达到预期范围或者迭代次数达到某个预设的阈值。在接收信号存在高斯白噪声时，可得三个LS拟合的代价函数分别为

其中，||·||_F表示矩阵的Frobenius范数。令及根据式(19)易知，对F₁、A₁和S₁的LS估计分别为

交替最小二乘算法在处理本发明所述三线性模型的具体步骤如下，a)初始化F₁、A₁和S₁；b)迭代计算和在计算过程中相应的F₁、A₁和S₁分别为上一次迭代后的结果；c)重复b直到迭代次数达到预设值或拟合误差达到预设阈值。

由于交替最小二乘算法在更新过程中F₁、A₁和S₁的误差将得到改善或者保持不变，但是不可能增大，因而交替最小二乘算法总是会收敛的。与大多数迭代算法，交替最小二乘算法所求得的极值仅是局部最小值。但是，在试验中所得到的极值就是全局最小值。此外，使用一些压缩算法可以进一步加快算法收敛。本发明在实际仿真中使用COMFAC算法，其主要是通过张量压缩的方法降低迭代计算的复杂度，一般仅需若干次迭代算法便可快速收敛。

定理1：对于式(17)中的PARAFAC模型，假设F₁、A₁和S₁的k-秩分别为k_F，k_A和k_S，若其满足

k_F+k_A+k_S≥2K+2 表达式(21)

则除了列模糊和尺度模糊之外，通过交替最小二乘获得的F₁、A₁和S₁的估计值是唯一的。列模糊和尺度模糊可以表示为

上式中Ω是一个列置换矩阵，N₁，N₂和N₃分别为误差矩阵，Δ₁，Δ₂和Δ₃为三个对角矩阵，且满足Δ₁Δ₂Δ₃＝I_K，I_K表示维数为K×K的矩阵。

F和A均为范德蒙矩阵，其具有旋转不变性。经过酉变换后的F₁和A₁仍然具有旋转不变特性，此外的旋转不变特性可表述为：

上式中，Ψ₁＝diag{-tan(πdsinθ₁/λ),，-tan(πdsinθ₂/λ),…,-tan(πdsinθ_K/λ)}，Ψ₂＝diag{-tan(πτ₁/L),-tan(πτ₂/L)…,-tan(πτ_K/L)}。J₁、J₂、J₃和J₄分别为：

其中，0表示元素全为0的矩阵，下标代表矩阵的维数。Re{·}、Im{·}分别为取实部和虚部。经过交替最小二乘，可获得F₁和A₁的估计值和根据式(23)的关系，可知-tan(πdsinθ_k/λ)与-tan(πτ_k/L)(k＝1,2,…,K)的LS估计值分别为

其中，与分别为和的第k列。最后，源信号的角度与时延可通过下式恢复

其中arcsin(·)与arctan(·)分别为反正弦与反正切运算。根据式(22)可知，与具有同步的列模糊特性，因此所估计的角度与时延是自动配对的。

对进行酉变换的运算量为2M²LDN+2M(LD)²N+8MLDN²。所提基于实数交替最小二乘的运算量为l[3O(K³)+6MLDNK+4MNK²+4NLDK²+MLDK²+2MNK+2LDNK+MLDK]/4，l为迭代的次数。恢复角度与时延的复杂度为[2K(M-1)(2M+1)+2K(LD-1)(LD+1)]/4。如表一所示，表一列举了无需谱峰搜索的ESPRIT与PARAFAC算法复杂度，可以看出，本发明算法复杂度低于PARAFAC，且在高维数据场景下(M或者LD很大时)，本发明算法和PARAFAC算法复杂度低于ESPRIT。

表一

由于可辨识度决定了算法最大可辨识参数数目，本发明算法的可辨识度由定理1决定。一般来说，k_A＝M，k_F＝LD，k_S＝2N，因而本发明算的最多可辨识(M+LD+2N)/2个信源，而PARAFAC算法最多可辨识(M+LD+N)/2个目标。因而在相同的条件下本发明算法可以识别更多的目标，且能够应对相干源场景。

进一步的，JADE的克拉美-罗界(Cram′er-Rao Bound，CRB)由下式给出

式中σ²为高斯白噪声的方差。U＝A⊙G，D＝[A′⊙G,A⊙G′]，A′和G′表示求导运算。

为验证本发明所述实值平行因子分解的多径参数估算方法的有效性，采用MATLAB数字仿真的形式对本发明所述实值平行因子分解的多径参数估算方法进行性能评估。仿真中假设K＝3多径，其角度与时延分别为(θ₁,τ₁)＝(10°,0.2s)(θ₂,τ₂)＝(20°,0.4s)、(θ₃,τ₃)＝(30°,0.6s)。接收天线数目为M＝6，L＝2，过采样率D＝4，数字序列数目为N＝16，多径信道衰落系数均设置为1。

图1是在SNR＝10dB非相干源条件下本发明算法进行200次蒙特卡洛仿真的散点图，图2是在SNR＝10dB、相干源(多径一和多径二的相干系数为1)条件下本发明算法200次蒙特卡洛实验的散点图。可以看出，两种仿真条件下三条多径的角度与时延参数均可被准确的估计出来。

为进一步分析本发明算法的估计性能，将ESPRIT算法、PARAFAC算法和本发明算法进行200次独立的蒙特卡洛仿真，参数估计的精度用均方根误差衡量(root mean squarederror，RMSE)，其具体定义为

式中，和r_k分别为在第i次仿真中估计的第条多径的角度(或时延)。

图3与图4分别为不同信噪比条件下本发明所述实值平行因子分解的多径参数估算方法对非相干源的角度和时延RMSE性能的对比。其中，图3为角度估计的RMSE曲线，图3为时延估计的RMSE曲线。由仿真结果可知，随着信噪比增加，所有算法的RMSE性能变好。本发明算法和PARAFAC算法性能优于ESPRIT，因为二者均利用了阵列数据的多维结构，而ESRPTI算法只利用了张量数据的一个维度的展开信息。此外，还可以看到本发明算法在低信噪比条件下性能优于PARAFAC算法，因为前后平滑技术或者了阵列数据的维数。尽管在高信噪比时PARAFAC算法更接近CRB，但注意到本发明算法的复杂度更低，因而本发明算法可在估计精度和计算复杂度方面达到良好的折衷。

图5与图6分别为所有算法在相干源场景下角度与时延估计的RMSE性能比较，其中第一个目标和第二个目标的相干度为1。此时，ESPRIT算法和PARAFAC算法均会失效，而本发明算法仍然能有效的工作。

本发明提出一种实值平行因子分解的多径参数估算方法。首先构建阵列数据的PARAFAC模型，然后利用前后平滑技术和酉变换技术对接收数据进行处理，构造阵列的增广输出的PARAFAC模型。再利用交替最小二乘算法获得相关导引矢量的估计，最后再利用阵列的旋转不变特性恢复源信号的角度与时延参数。本发明所提算法在进行交替最小二乘时仅涉及实数运算，故相比已有复数算法，本发明所提算法的计算复杂度更低，且由于前后平滑技术会扩展阵列接收数据，因而所提算法在低信噪比或者小快拍时能获得更优的估计性能，同时还能自动配对所估计角度和时延参数，并有效应对相干源场景，不需要对接收数据进行奇异值分解或者是谱峰搜索运算，在实际的无线通信场景中具有良好的应用前景。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种实值平行因子分解的多径参数估算方法，其特征在于，所述实值平行因子分解的多径参数估算方法包括如下步骤：

S1、构建接收数据的过采样矩阵，根据过采样矩阵构建出过采样矩阵阵列数据的PARAFAC模型；

S2、利用前后平滑技术和酉变换技术对过采样矩阵阵列数据的PARAFAC模型进行处理，构造过采样矩阵阵列的实数增广输出的PARAFAC模型；

S3、通过交替最小二乘算法获得实数增广输出的PARAFAC模型的导引矢量的估计；

S4、根据阵列的旋转不变特性，通过所述导引矢量的估计恢复源信号的角度与时延参数。

2.根据权利要求1所述实值平行因子分解的多径参数估算方法，其特征在于，所述步骤S2包括以下分步骤：

S21、利用前后平滑技术构建前后平滑后的PARAFAC模型张量；

S22、根据前后平滑后的PARAFAC模型张量构建扩展后的PARAFAC模型；

S23、利用酉变换技术对扩展后的PARAFAC模型进行实数变换接，构造过采样矩阵阵列的实数增广输出的PARAFAC模型。

3.根据权利要求2所述实值平行因子分解的多径参数估算方法，其特征在于，当利用M根均匀线性阵列的接收天线测量一个信道传输数字序列s_l，且接收信号x(t)以采样率D倍符号率被过采样，接收信号完全占据N个符号的时长，所述步骤S1中过采样矩阵如下：

式中T为符号周期；…,exp{-j2π(M-1)dsinθ_k/λ}]^T为接收导引矢量。

所述过采样矩阵阵列数据的PARAFAC模型如下：

式中，F即为时延矩阵；S为信源矩阵；表述维数为K×K×K的单位张量；A为方向矩阵。

4.根据权利要求3所述实值平行因子分解的多径参数估算方法，其特征在于，所述步骤S21中利用前后平滑技术构建前后平滑后的张量：

式中，Π_N表示反向交换矩阵，其下标表示矩阵的维数。

5.根据权利要求4所述实值平行因子分解的多径参数估算方法，其特征在于，所述步骤S22中扩展后的PARAFAC模型如下：

6.根据权利要求5所述实值平行因子分解的多径参数估算方法，其特征在于，所述步骤S23中通过酉变换将复数张量变换为实值张量如下：

式中， 为酉变换后的导引矢量，其下标表示矩阵的维数。

7.根据权利要求6所述实值平行因子分解的多径参数估算方法，其特征在于，所述步骤S3中相关导引矢量的估计如下：

式中，Z₁、Z₂和Z₃分别可被视为将张量数据沿着信源方向、时域方向和空域方向展开而获得的矩阵；及。

8.根据权利要求7所述实值平行因子分解的多径参数估算方法，其特征在于，所述步骤S4中所述源信号的角度与时延通过以下步骤恢复：

经过酉变换后的F₁和A₁仍然具有旋转不变特性，其旋转不变特性可表述如下：

其中，0表示元素全为0的矩阵，下标代表矩阵的维数；Re{·}、Im{·}分别为取实部和虚部；

经过交替最小二乘，获得F₁和A₁的估计值和如下：

式中，与分别为与的第k列；

则源信号的角度与时延可通过下式恢复：