CN103200136A - 一种频域超分辨率多径时延估计方法 - Google Patents

Abstract

本发明公开了一种频域超分辨率多径时延估计方法，包括如下步骤：（1）频域信道响应采样；（2）利用采样数据构造Hankle矩阵；（3）利用采样数据分块矩阵X₁和X₂构造矩阵E和矩E_I；（4）构造噪声子空间、模式矢量；（5）搜索多径时延。本发明提供的频域超分辨率多径时延估计方法，在相同信噪比的情况下，具有比传统的多重信号分类（MUSIC）算法还高的时间分辨率，并且不再需要奇异值分解或是特征值分解等计算量较大的矩阵运算操作，让多径时延估计更加实用、简单。

Description

一种频域超分辨率多径时延估计方法

技术领域

本发明涉及无线通信系统中的测距和定位技术，尤其涉及一种在信号频域进行超分辨率多径时延估计的方法。

背景技术

通常在无线传感器网络和蜂窝移动通信系统中，测距和定位技术广泛用于智能交通、生物医疗和工业控制等领域，实现目标定位、病情追踪和故障管理等应用。在所有的这些应用场景中，测距和定位的基础便是精确的时延估计技术，而在多径环境下多径时延的估计至关重要。

传统的超分辨率时延估计算法都是基于矩阵特征分解理论的，其本质是对相关矩阵特征分解后，依据信号子空间矩阵和噪声子空间矩阵的正交特性来得到多径信号的时延估计值。这些算法在信噪比较高的情况下时延估计分辨率都很高，但是，这些方法一般都需要奇异值分解或是特征值分解等计算量较大的矩阵运算操作。

发明内容

发明目的：为了克服现有技术中存在的不足，本发明提供一种频域超分辨率多径时延估计方法，能够在信噪比不高的情况下，使时延估计达到较高的分辨率，同时能够避免传统超分辨率时延估计算法的奇异值分解和特征值分解等运算，让多径时延估计更加实用、简单。

技术方案：为实现上述目的，本发明采用的技术方案为：

一种频域超分辨率多径时延估计方法，包括如下步骤：

第一步：频域信道响应采样

在频域利用矢量网络分析仪对信道频域响应进行N点的均匀采样，其采样频率范围为[f_L,f_H]，采样间隔记为△f=(f_H-f_L)/(N-1)，接收到的数据采样点分别为R(0),R(1),...,R(N-1)；其中f_L为频率下界，f_H为频率上界；

第二步：利用采样数据构造Hankle矩阵

将在频域采样得到的N个频域数据组织成为如下的L×(N-L+1)矩阵形式：（L为矩阵X的行数，比N小很多，并且可以变化，控制这计算复杂度和算法精确度）

其中，X₁和X₂是采样数据分块矩阵，且X₁是D×(N-L+1)矩阵，X₂是(L-D)×(N-L+1)矩阵，D是多径的传播路径数目，即：

第三步：利用采样数据分块矩阵X₁和X₂构造矩阵E和矩E_I

$E=\left[\begin{array}{c}{\left(X_1{X}_1^H\right)}^{-1}{X}_1{X}_2^H\\ -I\end{array}\right]$

$E_I=\left[\begin{array}{c}{\left({X_1}^H\right)}^{+}\\ -{\left({X_2}^{+}\right)}^H\end{array}\right]$

其中，矩阵I是(L-D)×(L-D)单位矩阵；

第四步：构造噪声子空间、模式矢量

构造噪声子空间矩阵Q=E×E^H或Q_I=E_I×E_I ^H以及模式矢量 $e=\left[\begin{array}{c}{e}^{-j2\mathrm{\π}{f}_L\mathrm{\τ}}{e}^{-j2\mathrm{\π}(f_L+\mathrm{\Δf})\mathrm{\τ}}...e^{-j2\mathrm{\π}(f_L+(L-1)\mathrm{\Δf})\mathrm{\τ}}\end{array}\right],$ 其中τ为传播路径的传播时延（j为序数符号，相当于虚数中的i）；

第五步：搜索多径时延

采用传统的MUSIC算法中的如下搜索方法来估计多径时延：

$\mathrm{MU}\left(\mathrm{\τ}\right)=\frac{1}{\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{\left|\right|e^H{q}_i\left|\right|}^2}$

其中，q_i是矩阵Q或者Q_I的第i列。

下面针对本案给出原理和演算说明。

假设信号r(t)是信号s(t)经过多个不同延时之后的叠加，如下式：

$r\left(t\right)=\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{\mathrm{\α}}_{d}s(t-{\mathrm{\τ}}_d)---\left(1\right)$

对上式进行傅里叶变换得到：

$R\left(f\right)={\∫}_{\∞}^{-\∞}r\left(t\right)e^{-j2\mathrm{\πft}}\mathrm{dt}$

$={\∫}_{-\∞}^{\∞}\left[\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{\mathrm{\α}}_{d}s(t-{\mathrm{\τ}}_d)\right]e^{-j2\mathrm{\πft}}\mathrm{dt}$

$=\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{\mathrm{\α}}_d{\∫}_{-\∞}^{\∞}s(t-{\mathrm{\τ}}_d)e^{-j2\mathrm{\πf}(t-{\mathrm{\τ}}_d)}{e}^{-j2\mathrm{\πf}{\mathrm{\τ}}_d}\mathrm{dt}---\left(2\right)$

$=\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{\mathrm{\α}}_d{e}^{-j2\mathrm{\πf}{\mathrm{\τ}}_d}{\∫}_{-\∞}^{\∞}s(t-{\mathrm{\τ}}_d)e^{-j2\mathrm{\πf}(t-{\mathrm{\τ}}_d)}\mathrm{dt}$

$=\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{\mathrm{\α}}_d{e}^{-j2\mathrm{\πf}{\mathrm{\τ}}_d}S\left(f\right)$

其中，S(f)是信号s(t)的傅里叶变换，即

在将式（2）的两边同时除以S(f)，得到：

$\frac{R\left(f\right)}{S\left(f\right)}=\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{\mathrm{\α}}_d{e}^{-j2\mathrm{\πf}{\mathrm{\τ}}_d}---\left(3\right)$

很明显，上式就是典型的频率估计问题。可以清楚地看到，其实时域上的时延估计问题经过傅里叶变换之后就会转换成频域上的频率估计问题。因此，很多超分辨率算法都被引进到时延估计问题中来。

在频域利用矢量网络分析仪可以建立如下多径传播信道模型：

$R\left(f\right)H\left(f\right)+W\left(f\right)=\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{\mathrm{\α}}_d{e}^{-j2\mathrm{\πf}{\mathrm{\τ}}_d}+W\left(f\right)---\left(4\right)$

其中，H(f)是相应频域中的多径信道响应函数，W(f)是零均值的加性高斯白噪声，D是多径的传播路径数目，α_d和τ_d分别是第d条路径的幅度和传播时延。

既然时域上的时延估计问题与频域上的频率估计问题等价，那么接下来的问题将会转换到频域上来说明。

不失一般性地，假设τ₁<τ₂...<τ_D，根据通常的多径传播情况，可以假设α₁>α₂...>α_D。另外，为了问题的简化，假设多径数目D是已知的。

对信道频域响应进行N点的均匀采样，其采样频率范围为[f_L,f_H]，采样间隔记为△f=(f_H-f_L)/(N-1)，接收到的数据采样点分别为R(0),R(1),...,R(N-1)；其中f_L为频率下界，f_H为频率上界；将式（4）离散化：

$R\left(k\right)=H\left(k\right)+W\left(k\right)=\underset{d=1}{\overset{D}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d{e}^{-j2\mathrm{\&pi;}{f}_k{\mathrm{\&tau;}}_d}+W\left(k\right)=\underset{d=1}{\overset{D}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_d{e}^{-j2\mathrm{\&pi;}{f}_L{\mathrm{\&tau;}}_d-j2\mathrm{\&pi;k\&Delta;f}{\mathrm{\&tau;}}_d}+W\left(k\right)---\left(5\right)$

PM-MUSIC算法

利用式（5）中接收到的数据采样点：R(0),R(1),...,R(N-1)，可以得到Hankle矩阵：

将式（6）改写成如下形式：

X=A(τ)×Φ(τ)×Κ+W        （7）

其中：

$K={\left[\begin{array}{cccc}{\mathrm{\&alpha;}}_1& {\mathrm{\&alpha;}}_2& \&CenterDot;\&CenterDot;\&CenterDot;& {\mathrm{\&alpha;}}_D\end{array}\right]}^T$

则：

其中：

U=[u₀,u₁,…,u_N-L]

u_i=[W(i),W(i+1),…,W(i+L-1)]^T

将式（8）中的A_L(τ)分成两个分块矩阵，即：

$A_L\left(\mathrm{\&omega;}\right)=\left[\begin{array}{c}{A}_{L1}\left(\mathrm{\&tau;}\right)\\ A_{L2}\left(\mathrm{\&tau;}\right)\end{array}\right]---\left(9\right)$

其中，A_L1(τ)是D×D矩阵，A_L2(τ)是(L-D)×D矩阵。

定义(L-D)×D传播算子矩阵P满足下式：

P^HA_L1(τ)=A_L2(τ)           （10）

同样地，将接收到的数据矩阵X和噪声矩阵U分别分成两个分块矩阵，即：

$X=\left[\begin{array}{c}{X}_1\\ X_2\end{array}\right]---\left(11\right)$

$U=\left[\begin{array}{c}{U}_1\\ U_2\end{array}\right]---\left(12\right)$

其中，X₁、U₁是D×(N-L+1)矩阵，X₂、U₂是(L-D)×(N-L+1)矩阵。

根据式（10）～式（12）可以得到：

将上式中的第一个等式两边同时乘以P^H，即：

根据式（13）和式（14）可以得到：

X₂-P^HX₁=U₂-P^HU₁         （15）

$\tilde{P}=\arg \min \left|\right|X_2-P^H{X}_1\left|\right|$        $\left(16\right)$

${=\left(X_1{X}_1^H\right)}^{-1}{X}_1{X}_2^H$

接下来，构造矩阵E如下：

$E=\left[\begin{array}{c}P\\ -I\end{array}\right]---\left(17\right)$

其中，矩阵I是(L-D)×(L-D)单位矩阵。

将上述构造的矩阵E的埃尔米特转置与A_L(τ)相乘，显然可以得到：

$E^H{A}_L\left(\mathrm{\&tau;}\right)=\left[\begin{array}{cc}{P}^H& -I\end{array}\right]\left[\begin{array}{c}{A}_{L1}\left(\mathrm{\&tau;}\right)\\ A_{L2}\left(\mathrm{\&tau;}\right)\end{array}\right]=P^H{A}_{L1}\left(\mathrm{\&tau;}\right)-A_{L2}\left(\mathrm{\&tau;}\right)=O---\left(18\right)$

在得到传播算子矩阵P估计值的前提下，构造的矩阵E的估计值便可得到，即：

$\tilde{E}=\left[\begin{array}{c}\tilde{P}\\ -I\end{array}\right]---\left(19\right)$

在加性高斯白噪声环境下，构造的矩阵E并不是正交的，于是引入正交投射矩阵Q来表示噪声子空间矩阵，即：

QA_L(τ)=O          （20）

其中：

$Q=\tilde{E}{\left({\tilde{E}}^H\tilde{E}\right)}^{-1}{\tilde{E}}^H---\left(21\right)$

在信号子空间矩阵A_L(ω)和噪声子空间矩阵Q都具备的情况下，利用MUSIC算法的搜索方法便可以得到下式：

$P_{\mathrm{MU}}\left(\mathrm{\&omega;}\right)=\frac{1}{\underset{i=0}{\overset{L-1}{\mathrm{\&Sigma;}}}{\left|\right|e^H{q}_i\left|\right|}^2}---\left(22\right)$

其中，q_i表示噪声子空间矩阵Q的第i列矢量，e表示模式矢量，即 $e=\left[\begin{array}{c}{e}^{-j2\mathrm{\&pi;}{f}_L\mathrm{\&tau;}}{e}^{-j2\mathrm{\&pi;}(f_L+\mathrm{\&Delta;f})\mathrm{\&tau;}}...e^{-j2\mathrm{\&pi;}(f_L+(L-1)\mathrm{\&Delta;f})\mathrm{\&tau;}}\end{array}\right].$

改进型PM-MUSIC算法

不同于PM-MUSIC算法的是，改进型PM-MUSIC算法构造新的矩阵E_I，并将其分解为两个分块矩阵E_I1和E_I2，即：

$E_I=\left[\begin{array}{c}{E}_{I1}\\ E_{I2}\end{array}\right]---\left(23\right)$

其中，E_I1和E_I2分别为D×(N-L+1)矩阵和(L-D)×(N-L+1)矩阵。

令E_I2是分块矩阵X₂广义逆的埃尔米特转置的负值，即E_I2=-(X₂ ⁺)^H，因此有：

E_I ^HA(τ)=E_I1 ^HA₁(τ)+E_I2 ^HA₂(τ)=E_I1 ^HA₁(τ)-X₂ ⁺A₂(τ)    （24）

其中，(·)⁺表示矩阵的广义逆。

从而矩阵E_I1可以被估计，即：

E_I1=argmin||E_I1 ^HX₁+E_I2 ^HX₂||²

=argmin||E_I1 ^HX₁-X₂ ⁺X₂||²     （25）

=(X₁X₁ ^H)⁺X₁(X₂ ⁺X₂)^H

=(X₁ ^H)⁺(X₁ ⁺X₁)(X₂ ⁺X₂)^H

根据矩阵广义逆的性质，当矩阵A∈C^m×n，则：

AA⁺=I_m的充分必要条件是rank(A)=m；

A⁺A=I_n的充分必要条件是rank(A)=n。

所以，当L≥(N+D+1)/2时，rank(X₂)=N-L+1，可以得到X₂ ⁺X₂=I，否则X₂ ⁺X₂≠I；同样地，当L≥N-D时，rank(X₁)=N-L+1，可以得到X₁ ⁺X₁=I，否则X₁ ⁺X₁≠I。在不过多降低性能的情况下，为了进一步降低计算复杂度，可以令E_I1=(X₁ ^H)⁺。

这样一来就获得了矩阵E_I的估计值，再构造噪声子空间矩阵Q_I=E_I×E_I ^H和模式矢量 $e=\left[\begin{array}{c}{e}^{-j2\mathrm{\&pi;}{f}_L\mathrm{\&tau;}}{e}^{-j2\mathrm{\&pi;}(f_L+\mathrm{\&Delta;f})\mathrm{\&tau;}}...e^{-j2\mathrm{\&pi;}(f_L+(L-1)\mathrm{\&Delta;f})\mathrm{\&tau;}}\end{array}\right].$

最后，同样采用传统的MUSIC算法中的如下搜索方法来估计多径时延：

$\mathrm{MU}\left(\mathrm{\&tau;}\right)=\frac{1}{\underset{i=0}{\overset{L-1}{\mathrm{\&Sigma;}}}{\left|\right|e^H{q}_i\left|\right|}^2}---\left(26\right)$

其中，q_i是矩阵Q或者Q_I的第i列。

有益效果：本发明提供的频域超分辨率多径时延估计方法，在相同信噪比的情况下，具有比传统的多重信号分类（MUSIC）算法还高的时间分辨率，并且不再需要奇异值分解或是特征值分解等计算量较大的矩阵运算操作，让多径时延估计更加实用、简单。

附图说明

图1为城市多径效应示意图；

图2为本发明方法的流程图；

图3a为改进型PM-MUSIC与PM-MUSIC和MUSIC多径时延估计对比图：τ₁=20ns,τ₂=40ns,τ₃=60ns,τ₄=80ns；

图3b为改进型PM-MUSIC与PM-MUSIC和MUSIC多径时延估计对比图：τ₁=20ns,τ₂=70ns,τ₃=120ns,τ₄=170ns；

图4a为改进型PM-MUSIC与PM-MUSIC的MSE对比图：信噪比；

图4b为改进型PM-MUSIC与PM-MUSIC的MSE对比图：接收数据数目N；

图4c为改进型PM-MUSIC与PM-MUSIC的MSE对比图：接收数据Hankle矩阵行数L。

具体实施方式

下面结合附图对本发明作更进一步的说明。

如图1所示为城市多径效应的示意图，针对这种多径效应的一种频域超分辨率多径时延估计方法流程图如图2所示，包括如下步骤：

第一步：频域信道响应采样

在频域利用矢量网络分析仪对信道频域响应进行N点的均匀采样，其采样频率范围为[f_L,f_H]，采样间隔记为△f=(f_H-f_L)/(N-1)，接收到的数据采样点分别为R(0),R(1),...,R(N-1)；其中f_L为频率下界，f_H为频率上界；

第二步：利用采样数据构造Hankle矩阵

将在频域采样得到的N个频域数据组织成为如下的L×(N-L+1)矩阵形式：

其中，X₁和X₂是采样数据分块矩阵，且X₁是D×(N-L+1)矩阵，X₂是(L-D)×(N-L+1)矩阵，D是多径的传播路径数目，即：

（28）

第三步：利用采样数据分块矩阵X₁和X₂构造矩阵E和矩E_I

$E=\left[\begin{array}{c}{\left(X_1{X}_1^H\right)}^{-1}{X}_1{X}_2^H\\ -I\end{array}\right]$       （29）

$E_I=\left[\begin{array}{c}{\left({X_1}^H\right)}^{+}\\ -{\left({X_2}^{+}\right)}^H\end{array}\right]$

其中，矩阵I是(L-D)×(L-D)单位矩阵；

第四步：构造噪声子空间、模式矢量

构造噪声子空间矩阵Q=E×E^H或Q_I=E_I×E_I ^H以及模式矢量 $e=\left[\begin{array}{c}{e}^{-j2\mathrm{\&pi;}{f}_L\mathrm{\&tau;}}{e}^{-j2\mathrm{\&pi;}(f_L+\mathrm{\&Delta;f})\mathrm{\&tau;}}...e^{-j2\mathrm{\&pi;}(f_L+(L-1)\mathrm{\&Delta;f})\mathrm{\&tau;}}\end{array}\right]$ 其中τ为传播路径的传播时延；

第五步：搜索多径时延

采用传统的MUSIC算法中的如下搜索方法来估计多径时延：

$\mathrm{MU}\left(\mathrm{\&tau;}\right)=\frac{1}{\underset{i=0}{\overset{L-1}{\mathrm{\&Sigma;}}}{\left|\right|e^H{q}_i\left|\right|}^2}---\left(30\right)$

其中，q_i是矩阵Q或者Q_I的第i列。

仿真结果分析

在以下的仿真实验中，考虑有四条多径，其幅度分别为α₁=1,α₂=0.8,α₃=0.6,α₄=0.4，为了进一步说明传统MUSIC算法、PM-MUSIC算法和改进型PM-MUSIC算法的时间分辨率，在不同实验中的时延是并不相同的，另外，频域采样的数据中频率下界f_L为1GHz，采样频率间隔△f为2MHz。

如图3（图3a和图3b）所示，仿真实验比较改进型PM-MUSIC算法与传统MUSIC算法以及PM-MUSIC算法的时间分辨率，实验是在高斯白噪声环境下信噪比为20dB的情况下进行的。接收数据数目N为200，传统MUSIC算法中相关矩阵的维数P取20，而PM-MUSIC算法和本案提出的改进型PM-MUSIC算法中的L为25。当τ₁=20ns,τ₂=40ns,τ₃=60ns,τ₄=80ns时，传统MUSIC算法有的尖峰不在正确的时延位置上，而且有的位置上的尖峰非常小；PM-MUSIC算法除了在第一个时延点上有明显的尖峰外，其余的时延点上没有尖峰；但是，本案提出的改进型PM-MUSIC算法在每个时延点上都有明显的尖峰。也就是说，前两者的时间分辨率比较低，不能很好的区分各条多径时延，而改进型PM-MUSIC算法能很好的区分每条多径时延。当τ₁=20ns,τ₂=70ns,τ₃=120ns,τ₄=170ns时，虽然三种算法在多径时延点上都有尖峰，但是不难发现，除了在第一个时延点上三种算法的尖峰类似以外，在其他的时延点上都是PM-MUSIC算法的尖峰优于传统MUSIC算法，而改进型PM-MUSIC最优。

如图4（图4a、图4b和图4c）所示，分别进行了1000次独立的Monte-Carlo仿真实验，比较PM-MUSIC算法与改进型PM-MUSIC算法在不同参数下的性能。定义MSE为：

${\mathrm{MSE}}_{\mathrm{dB}}=10{\log }_{10}\left(\frac{1}{\mathrm{ND}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{D}{\mathrm{\&Sigma;}}}{({\mathrm{\&tau;}}_{\mathrm{ij}}-{\widetilde{\mathrm{\&tau;}}}_{\mathrm{ij}})}^2\right)---\left(31\right)$

不难发现，当信噪比高于0dB时，改进型PM-MUSIC算法要优于PM-MUSIC算法。另外，随着接收数据数目N的不断加大，PM-MUSIC算法的性能提升有限，而改进型的PM-MUSIC提升的性能更为明显，当N等于400的时候该算法的均方误差要小8dB左右。随着接收数据矩阵行数L的不断增加，PM-MUSIC算法与改进型的PM-MUSIC的性能都会有一定程度的改善，但是PM-MUSIC的性能提升得更多。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (1)

1.一种频域超分辨率多径时延估计方法，其特征在于：包括如下步骤：

第一步：频域信道响应采样

在频域利用矢量网络分析仪对信道频域响应进行N点的均匀采样，其采样频率范围为[f_L,f_H]，采样间隔记为△f=(f_H-f_L)/(N-1)，接收到的数据采样点分别为R(0),R(1),...,R(N-1)；其中f_L为频率下界，f_H为频率上界；

第二步：利用采样数据构造Hankle矩阵

将在频域采样得到的N个频域数据组织成为如下的L×(N-L+1)矩阵形式：

其中，X₁和X₂是采样数据分块矩阵，且X₁是D×(N-L+1)矩阵，X₂是(L-D)×(N-L+1)矩阵，D是多径的传播路径数目，即：

第三步：利用采样数据分块矩阵X₁和X₂构造矩阵E和矩E_I

$E=\left[\begin{array}{c}{\left(X_1{X}_1^H\right)}^{-1}{X}_1{X}_2^H\\ -I\end{array}\right]$

$E_I=\left[\begin{array}{c}{\left({X_1}^H\right)}^{+}\\ -{\left({X_2}^{+}\right)}^H\end{array}\right]$

其中，矩阵I是(L-D)×(L-D)单位矩阵；

第四步：构造噪声子空间、模式矢量

构造噪声子空间矩阵Q=E×E^H或Q_I=E_I×E_I ^H以及模式矢量 $e=\left[\begin{array}{c}{e}^{-j2\mathrm{\&pi;}{f}_L\mathrm{\&tau;}}{e}^{-j2\mathrm{\&pi;}(f_L+\mathrm{\&Delta;f})\mathrm{\&tau;}}...e^{-j2\mathrm{\&pi;}(f_L+(L-1)\mathrm{\&Delta;f})\mathrm{\&tau;}}\end{array}\right],$ 其中τ为传播路径的传播时延；

第五步：搜索多径时延

采用传统的MUSIC算法中的如下搜索方法来估计多径时延：

$\mathrm{MU}\left(\mathrm{\&tau;}\right)=\frac{1}{\underset{i=0}{\overset{L-1}{\mathrm{\&Sigma;}}}{\left|\right|e^H{q}_i\left|\right|}^2}$

其中，q_i是矩阵Q或者Q_I的第i列。

Images