CN102208930B - 快速信号子空间估计方法 - Google Patents

Abstract

本发明涉及数字信号处理技术领域。本发明提供了一种快速信号子空间估计方法，所述方法包括：根据观测数据模型，获得接收天线阵元的初始观测信号；构造新的观测数据，对所述新的观测数据进行迭代；根据所述迭代结果构造转换矩阵；以及，获得所述转换矩阵的列向量。本发明的方法在数据层次上进行多级维纳滤波前向迭代，得到前向迭代过程的转换矩阵，进而利用Householder矩阵的酉性质快速地计算该转换矩阵，从而获得信号子空间的一组标准基向量；所述方法借助了多级维纳滤波前向迭代的Householder变换，不仅进一步降低了计算复杂度，还保证了算法在有限精度下的数值稳健性。

Description

快速信号子空间估计方法

技术领域

本发明涉及数字信号处理技术领域，特别涉及一种降低信号子空间估计计算复杂度并保证数值稳健性的方法。

背景技术

随着大规模集成电路和数字信号处理技术的发展，信号子空间估计问题已经成为许多领域的关键问题之一，如相控阵雷达的波达方向(DOA)估计和降秩空时自适应处理(STAP)、移动通信的多用户检测(Multi-user Detection)等。传统的信号子空间估计方法通常是对观测数据协方差矩阵进行特征值分解(EVD)，对应大特征值的特征向量张成信号子空间，其余特征向量张成噪声子空间。然而由于特征值分解的计算复杂度非常高，数量级为O(M³)，其中M为观测数据的维数，并且需要计算观测数据协方差矩阵，计算复杂度为O(M²N)，其中N为形成协方差矩阵的采样支持长度。因此利用特征值分解估计信号子空间的计算复杂度为O(M³+M²N)。然而在很多实际应用中，特别是相控阵雷达和第三代移动通信系统，观测数据维数M通常很大，因此巨大的计算量往往限制了实时应用。

学术界和产业界相继提出了一系列的信号子空间快速估计方法来避免计算特征值分解，如G.Xu提出了基于Lanczos的快速子空间估计方法，其计算复杂度为O(M²P+M²N)，其中P为信号子空间的维数。当P远小于M时，该方法的计算量得到一定程度减轻，然而其仍然需要计算协方差矩阵。另外一种方法是基于多级维纳滤波(MSWF)前向迭代来估计信号子空间。多级维纳滤波前向迭代可以在数据层次(data-level)或者协方差矩阵层次(power-level)实现。协方差矩阵层次的前向迭代与G.Xu的方法的计算量一致，而数据层次的前向迭代，特别是基于相关相减结构的阻塞矩阵选择，能够大大降低其计算复杂度，然而其数值稳健性无法得到保证，尤其在有限精度运算中。

发明内容

本发明的目的是克服上述缺陷，提出一种数值稳健的快速信号子空间估计方法。

为实现上述目的，本发明的快速子空间估计方法包括：

获得天线阵元的观测数据模型；

构造新的观测数据，对所述新的观测数据进行迭代；

根据所述迭代结果构造转换矩阵；以及

获得所述转换矩阵的列向量。

进一步地，所述方法在所述构造新的观测数据，对所述新的观测数据进行迭代的步骤之前还包括：

根据所述观测数据模型获得天线阵元接收信号数据；

根据所述天线阵元接收信号数据构造新的期望信号和新的观测数据；

其中，所述观测数据模型为：

x(k)＝[x₁(k)，x₂(k)，…，x_M(k)]^T，k＝1，2，…，N

＝A(θ)s(k)+n(k)

＝[a(θ₁)，a(θ₂)，…，a(θ_P)][s₁(k)，s₂(k)，…，s_P(k)]^T

+[n₁(k)，n₂(k)，…，n_M(k)]^T

其中，x_M(k)代表第M个天线阵元的信号数据；s(k)，A(θ)和n(k)分别表示信号矢量、导向矩阵以及噪声矢量；N为采样支持长度，a(θ_i)代表波达方向为θ_i的信号导向矢量，并且其具有如下形式：

其中

d和λ分别表示ULA阵的阵元间距以及接收信号的载波波长；

进一步地，所述对所述新的观测数据进行迭代的步骤通过多级维纳滤波前向迭代计算实现；

进一步地，所述迭代结果为：归一化互相关矢量h_i和阻塞矩阵Q_B，i；

进一步地，所述构造后的转换矩阵为：

$T_P=[{\mathrm{\κ}}_1{q}_1,{\mathrm{\κ}}_2{Q}_{B,1}^H{q}_2,\·\·\·,{\mathrm{\κ}}_P\left(\underset{k=1}{\overset{P-1}{\mathrm{\Π}}}{Q}_{B,k}^H\right)q_P]$

$=[h_1,Q_{B,1}^H{h}_2,\·\·\·,\left(\underset{k=1}{\overset{P-1}{\mathrm{\Π}}}{Q}_{B,k}^H\right)h_P]$

其中，k_i为具有单位范数，任意相位的旋转因子；h_i为归一化互相关矢量；Q_Bi为阻塞矩阵；q_i为Householder矩阵Q_i的第一列；

进一步地，所述获得所述转换矩阵的列向量的步骤通过按照多级维纳滤波前向迭代的Householder矩阵后向连乘来实现；

进一步地，所述新的期望信号为天线参考阵元的观测数据，即d₀(k)＝x₁(k)；所述新的观测数据为：x₀(k)＝[x₂(k)，x₃(k)，…，x_M(k)]^T，其中，x_M(k)代表第M个天线阵元的信号数据；

进一步地，所述多级维纳滤波前向迭代计算在所述新的观测数据的数据层次上进行。

本发明的方法不需要计算观测数据协方差矩阵或者特征值分解，而是在数据层次上进行多级维纳滤波前向迭代，得到前向迭代过程的转换矩阵，进而利用Householder矩阵的酉性质快速地计算该转换矩阵，从而获得信号子空间的一组标准基向量；所述方法借助了多级维纳滤波前向迭代的Householder变换，不仅进一步降低了计算复杂度，还保证了算法在有限精度下的数值稳健性。

附图说明

图1是本发明具体实施例的ULA天线阵观测数据模型示意图；

图2是本发明具体实施例的基于Householder变换的多级维纳滤波前向迭代的示意图；

图3是M＝20，K＝19，N＝4K时本发明具体实施例与传统算法的乘法计算量比较示意图；

图4是M＝40，K＝39，N＝4K时本发明具体实施例与传统算法的乘法计算量比较示意图；

具体实施方式

本发明提出的快速子空间估计方法，结合附图和实施例说明如下。

S101：获得天线阵元的观测数据模型；

M个单元的均匀天线阵(ULA)接收到P个窄带远场信号，在第k个快拍时刻，M×1维的复观测数据可以表示为：

x(k)＝[x₁(k)，x₂(k)，…，x_M(k)]^T，k＝1，2，…，N

＝A(θ)s(k)+n(k)    (1)

＝[a(θ₁)，a(θ₂)，…，a(θ_P)][s₁(k)，s₂(k)，…，s_P(k)]^T

+[n₁(k)，n₂(k)，…，n_M(k)]^T

其中，x_M(k)代表第M个天线阵元的信号数据；s(k)，A(θ)和n(k)分别表示信号矢量、导向矩阵以及噪声矢量；N为采样支持长度，a(θ_i)代表波达方向为θ_i的信号导向矢量，并且其具有如下形式：

其中

d和λ分别表示ULA阵的阵元间距以及接收信号的载波波长。

首先，对于接收到的信号，假设其为窄带信号，具有不同的波达方向θ_i，并且所有信号为零均值的联合平稳高斯过程，即E{s(k)}＝0；此外，对于噪声，假设其为时域白、零均值的平稳高斯过程，即E{n(k)}＝0；并且有：信号之间不相关，信号和噪声之间不相关，不同阵元的噪声之间亦不相关，即有E{s(k)s^H(k)}＝R_s，E{s(k)n^H(k)}＝0，E{n(k)n^H(k)}＝σ²I。由此，式(1)构成了天线阵元接收信号的基本观测模型。

如图1所示为本发明实施例的天线阵接收数据的观测数据模型；

S102：构造新的观测数据，对所述新的观测数据进行迭代；

更具体地：

S1021：根据所述观测数据模型获得天线阵元接收信号数据；

S1022：根据所述天线阵元接收信号数据构造新的期望信号和新的观测数据；

S1023：对所述新的观测数据进行迭代。

本实施例采用基于Householder变换的多级维纳滤波前向迭代，所述迭代过程基于数据层次实现，避免了观测数据协方差矩阵的计算，从而降低了计算复杂度。

在具体实施过程中，为了能够执行多级维纳滤波前向迭代，在进行迭代计算之前，必须通过各个阵元接收的观测数据得到先验的本地期望信号(也称为训练信号)d₀(k)以及输入数据向量x₀(k)；在本实施例中，将d₀(k)设置为天线参考阵元的观测数据x₁(k)(第一个阵元数据)，即d₀(k)＝x₁(k)；而多级维纳滤波前向迭代的输入数据向量x₀(k)则表示为x₀(k)＝[x₂(k)，x₃(k)，…，x_M(k)]^T，迭代计算在此新的观测数据基础上进行。

经过多次的前向迭代后，可以得到每级互相关矢量r_i、范数归一化互相关矢量h_i、Householder矩阵Q_i、Householder向量q_i和阻塞矩阵Q_Bi，它们之间存在如下的数学关系，其中k_i是具有单位范数的复相位旋转因子，即，即||k_i||₂＝1：

$Q_i=[q_i,Q_{B,i}^H],q_i={\mathrm{\κ}}_i^{*}{h}_i$

h_i＝r_i/||r_i||₂    (3)

Q_B，ih_i＝Q_B，iq_i＝0

在基于Householder变换的多级维纳滤波前向迭代中，利用h_i、||k_i||₂＝1以及

的关系，可计算出各级的阻塞矩阵Q_B，i。

如图2所示，在具体实施过程中，可以通过以下方法完成基于Householder变换的多级维纳滤波前向迭代：

for i＝1:P

$r_i=E\left\{x_{i-1}\left(k\right)d_{i-1}^{*}\left(k\right)\right\}$

${\mathrm{\δ}}_i=\sqrt{r_i^H{r}_i}$

h_i＝r_i/δ_i

k_i＝-sign(real(h_i，1))

v_i＝h_i-k_iδ_iu_i    (4)

${\mathrm{\β}}_i=-1/\left({\mathrm{\κ}}_i^{*}{\mathrm{\δ}}_i{v}_{i,1}\right)$

end

从i＝1到P执行上述过程，共进行P次迭代。

在迭代过程中，Householder矩阵Q_i不需要显式地计算出来，即无需计算

S103：根据所述迭代的结果构造转换矩阵；

假设多级维纳滤波前向迭代进行P次后得到的转换矩阵由T_P表示，它直接作用于输入的期望信号d₀(k)和输入数据向量x₀(k)，从而可得到一组期望信号向量d(k)＝[d₀(k)，d₁(k)，…，d_P-1(k)]^T，并且期望信号向量d(k)的协方差矩阵为三对角矩阵，即

其中R_x和R_d分别表示输入数据向量x₀(k)和输出期望信号向量d(k)的协方差矩阵，且R_d为三对角矩阵。基于步骤S102计算得到的归一化互相关矢量h_i和阻塞矩阵Q_B，i、图2所示的多级维纳滤波前向迭代过程以及前述的||k_i||₂＝1，转换矩阵T_P可以表示如下：

$T_P=[{\mathrm{\κ}}_1{q}_1,{\mathrm{\κ}}_2{Q}_{B,1}^H{q}_2,\·\·\·,{\mathrm{\κ}}_P\left(\underset{k=1}{\overset{P-1}{\mathrm{\Π}}}{Q}_{B,k}^H\right)q_P]$ (5)

$=[h_1,Q_{B,1}^H{h}_2,\·\·\·,\left(\underset{k=1}{\overset{P-1}{\mathrm{\Π}}}{Q}_{B,k}^H\right)h_P]$

其中，其中，k_i为具有单位范数，任意相位的旋转因子；h_i为归一化互相关矢量；Q_B，i为阻塞矩阵；q_i为Householder矩阵Q_i的第一列。

T_P的列向量张成了信号子空间，并且这组列向量是信号子空间的一组标准正交基，所以，通过计算所述转换矩阵并获得其列向量即可实现信号子空间估计。

在基于Householder变换的多级维纳滤波前向迭代中，Householder矩阵Q_i具有严格的酉性质，即

因此酉变换减小了有限精度运算下的舍入误差在前向迭代过程中的传播，T_P列向量(即信号子空间的标准正交基)的正交性得到了非常好的保证，理论上可以证明其正交性误差非常接近机器的运算精度，即有限精度运算下的正交性误差与2^-ε在同一个数量级，其中ε为定点数中表示小数部分的比特长度，因此子空间估计的数值稳健性得到显著增强。

S104：获得所述转换矩阵的列向量；

本实施例采用转换矩阵T_P的快速计算方法获得其列向量，从而大大减小了信号子空间估计的计算量。利用Householder矩阵Q_i的酉性质

从数学上可以证明：转换矩阵T_P是矩阵Q的前P列，并且Q可表示为如下一系列矩阵

的连乘：

$Q=[q_1,Q_{B,1}^H]\left[\begin{array}{cc}1& 0\\ 0& [q_2,Q_{B,2}^H]\end{array}\right]\·\·\·\left[\begin{array}{cc}{I}_{P-2}& 0\\ 0& [q_{P-1},Q_{B,P-1}^H]\end{array}\right]$

$=\underset{k=1}{\overset{P-1}{\mathrm{\Π}}}\left[\begin{array}{cc}{I}_{k-1}& 0\\ 0& Q_k\end{array}\right]---\left(6\right)$

$=\underset{k=1}{\overset{P-1}{\mathrm{\Π}}}{\stackrel{\‾}{Q}}_k$

其中矩阵

的(k-1)×(k-1)主子式为单位矩阵I_k-1。计算矩阵

有两种方法：前向连乘和后向连乘。在矩阵

的后向连乘中(即第一步先计算

，再计算

，依此类推)，中间矩阵的(k-1)×(k-1)主子式逐渐消失。利用该性质，本实施例的转换矩阵T_P的快速计算方法为：不按照多级维纳滤波前向迭代的顺序前向连乘来计算矩阵

(即第一步先计算，再计算

，依此类推)，而是通过其后向连乘来进一步降低计算量。

具体实施过程中，可以通过以下方法完成转换矩阵T_P的快速计算：

初始化：Q＝I

连乘：for i＝(P-1):1

$Q=Q_i(i:K,i:K)\×Q(i:K,i:K)$                                   (7)

$=(I-{\mathrm{\β}}_i{v}_i(i:K)v_i^H(i:K))\×Q(i:K,i:K)$

end

截取：T_P＝Q(:，1:P)

其中，K＝M-1，I为单位矩阵；Q_i为MSWF前向迭代过程中的Householder矩阵；v_i为Householder向量；β_i为系数因子。

通过上述方法得到转换矩阵T_P的计算结果，从而获得其列向量，所述列向量即为估计信号子空间的一组标准正交基。

传统的基于相关相减结构的多级维纳滤波亦可得到信号子空间的一组标准正交基，并且这组基向量即为范数归一化互相关矢量h_i。然而在基于相关相减结构的多级维纳滤波前向迭代中，阻塞矩阵B_i不具有酉性质，即

将导致有限精度运算下的舍入误差通过多级维纳滤波前向迭代而传播，舍入误差产生积累，最终使得信号子空间基向量的正交性弱化，算法产生数值不稳健。本发明给出的方法中，数值稳健性通过Householder变换得到非常好的保证，因为Householder变换具有酉变换性质，多级维纳滤波的前向迭代不传播舍入误差，并且信号子空间基向量的正交性得到保证，非常接近于机器的运算精度。此外，通过步骤S104的转换矩阵T_P的快速计算方法，其计算量亦可以得到进一步降低。

表1所示为本发明的计算方法、传统相关相减结构的多级维纳滤波方法以及基于公式(5)直接计算转换矩阵T_P这三种方法的乘法计算量比较(实际应用中通常不考虑加法计算量)，其中K＝M-1。表1的乘法计算量公式难以明显比较三种方法的计算量。

表1

如图3、4所示为具体实施过程中，当M＝20或40(K＝19或39)以及N＝4K时本发明的方法与其余两种传统方法的乘法计算量比较。由图3和4可以看出，本发明提出的信号子空间估计方法，可以大大减小计算量，尤其对于天线阵接收信号较多、多级维纳滤波前向迭代次数大的环境。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (7)

1.一种快速信号子空间估计方法，其特征在于，包括： 

获得天线阵元的观测数据模型；其中，所述观测数据模型为： 

x(k)=[x₁(k),x₂(k),…，x_M(k)]^T,k=1,2,…,N 

=A(θ)s(k)+n(k) 

=[a(θ₁),a(θ₂),…,a(θ_P)][s₁(k),s₂(k),…,s_P(k)]^T

+[n₁(k),n₂(k),…,n_M(k)]^T

其中，x_M(k)代表第M个天线阵元的信号数据；s(k)，A(θ)和n(k)分别表示信号矢量、导向矩阵以及噪声矢量；N为采样支持长度，a(θ_i)代表波达方向为θ_i的信号导向矢量，并且其具有如下形式： 

其中

d和λ分别表示ULA阵的阵元间距以及接收信号的载波波长;

构造新的观测数据，对所述新的观测数据进行迭代；其中，对所述新的观测数据进行迭代的步骤通过多级维纳滤波前向迭代计算实现； 

根据所述迭代结果构造转换矩阵；以及 

获得所述转换矩阵的列向量；其中，构造后的转换矩阵为： 

其中，κ_i为具有单位范数，任意相位的旋转因子；h_i为归一化互相关矢量；Q_B,i为阻塞矩阵；q_i为Householder矩阵Q_i的第一列。 

2.如权利要求1所述的方法，其特征在于，所述方法在构造新的观测数据、对所述新的观测数据进行迭代的步骤之前还包括： 

根据所述观测数据模型获得天线阵元的初始接收信号数据； 

根据所述天线阵元的初始接收信号数据构造新的期望信号和新 的观测数据。 

3.如权利要求1所述的方法，其特征在于，所述迭代结果为：归一化互相关矢量h_i和阻塞矩阵Q_B,i。 

4.如权利要求1所述的方法，其特征在于，所述获得所述转换矩阵的列向量的步骤通过按照多级维纳滤波前向迭代的Householder矩阵后向连乘来实现。 

5.如权利要求2所述的方法，其特征在于，所述新的期望信号为天线参考阵元的观测数据，即d₀(k)=x₁(k)。 

6.如权利要求2所述的方法，其特征在于，所述新的观测数据为：x₀(k)=[x₂(k),x₃(k),…,x_M(k)]^T，其中，x_M(k)代表第M个天线阵元的信号数据。 

7.如权利要求1所述的方法，其特征在于，所述多级维纳滤波前向迭代计算在所述新的观测数据的数据层次上进行。 

Images