CN102662158B - 一种对传感器天线阵列接收信号的快速处理方法 - Google Patents

Abstract

该发明属于对传感器天线阵列接收信号的快速处理方法。包括：划分子阵及建立各子阵的信号矩阵，建立自相关矩阵和互相关矩阵，建立信号子空间矩阵，建立信号特征矩阵及其逆矩阵，建立过渡矩阵，建立旋转不变关系矩阵，最后获得不同信号的来波方向角。该发明采用直接建立旋转不变关系矩阵的方法，只利用矩阵的乘法和加法处理构造阵列旋转不变关系的方法实现传感器天线阵列对单个或多个无线传播信号的来波方向进行快速测定；与背景技术相比不但处理速度快，特别是在信噪比较低时对来波方向角的测定精度显提高，而具有对设备和硬件性能要求低，对目标来波方向的测定快速、准确，精度高，流程短、运行成本低，应用范围广等特点。

Description

一种对传感器天线阵列接收信号的快速处理方法

技术领域

本发明属于利用传感器天线阵列对无线传播信号的来波方向进行测定的领域，特别是一种利用天线阵列的旋转不变性（ESPRIT）对单个或多个无线传播信号的来波方向进行快速测定（定向）的方法。采用该方法可在无需为了得到旋转不变关系矩阵进行多次特征值分解的情况下进行快速测定（定向）、并可获得较好的效果（精确度），以及降低测定成本等优点。

背景技术

接收无线传播信号的天线阵列信号处理技术已广泛应用。阵列信号处理是将一组传感器在空间的不同位置按照一定的规则布置，形成传感器阵列，通过多个通道接收，利用一定的空域处理或者空时等多维联合处理的方式，得到所需的信号或信号参数。与传统的单个定向传感器相比，阵列信号处理具有更灵活的波束控制，更高的信号增益和空间超分辨能力，与此相关的研究工作迅速发展，应用范围也不断扩大，现已成功应用于地震、天文、雷达、声纳、医学信号检测、无线通信、地质勘探等领域。信号的来波方向(DOA)是阵列信号处理技术所需要的一个十分重要的参数，是阵列信号处理领域研究的热点方向。现代超分辨的信号来波方向估计（测定）的基本原理是通过将阵列多个阵元在同一时刻的采样信号排列为一向量，计算多个快拍下的阵列接收信号的自相关矩阵，并建立相应的信号子空间和噪声子空间；利用信号子空间的特性进行来波方向角的估计（测定）。多重信号分类法(MUSIC算法)是超分辨测向技术的代表之一，它利用信号子空间和噪声子空间的正交性原理构造一个伪谱函数，并在一定的方位内进行谱峰搜索，虽然多重信号分类法能得到很高的角度分辨率，但是其谱峰搜索尤其是在高精度的二维搜索时需要的庞大的计算量阻碍了其应用。而以旋转不变技术（ESPRIT）为代表的方法利用局部阵列的平移不变性，无需计算量庞大的谱峰搜索就能得到信号来波方向角的超分辨估计（测定），其估计（测定）的精度与MUSIC（多重信号分类法）相差无几。因此利用阵列的旋转不变性来估计信号的来波方向角具有重大的研究意义。在ESPRIT处理方法的具体实现方面，典型的有基于不同准则的ESPRIT算法和基于ESPRIT的推广算法。其中，基于最小二乘准则的LS-ESPRIT算法和基于总体最小二乘准则的TLS-ESPRIT法是两种基于不同准则的ESPRIT算法。LS-ESPRIT算法是用一个范数平方为最小的扰动去干扰信号子空间，校正信号子空间中的噪声的方法。TLS-ESPRIT算法是同时扰动两个子阵的信号子空间，并使扰动范数的平方保持最小，同时校正两个子阵的信号子空间中存在的噪声。这两种算法都需要对由两个子阵接收信号组成的组信号向量的协方差矩阵进行特征值分解，亦即要对2倍于子阵的协方差矩阵阶数的矩阵进行特征值分解，因而计算量庞大；而且，在二者构造旋转不变关系矩阵的过程中都需要再次对计算过程中构造的矩阵进行特征值分解；而后者在计算中还需要对一般矩阵进行求逆的计算，这样导致了ESPRIT算法实施需要的计算量与需要谱峰搜索的MUSIC算法比起来也相差无几；而矩阵束ESPRIT算法则是利用两个子阵接收数据的自相关矩阵和互相关矩阵进行广义特征分解以确定来波方向的方法，该方法虽然可以有效的降低ESPRIT算法的运算量，但是对来波方向测量精度难以满足常规应用的要求。

发明内容

本发明的目的是针对背景技术存在的缺陷，研究设计一种对传感器天线阵列接收信号的快速处理方法，以克服背景技术存在的缺陷，达到在降低对设备和硬件性能要求的同时，对目标来波方向进行快速、准确测定，提高精度，简化流程、降低运行成本等目的。

本发明的解决方案是首先将均匀直线阵列按照图1所示的划分方法分为两个子阵，将第一个子阵得到的采样信号矩阵与自身进行相关运算得到子阵1的自相关矩阵，再将第一个子阵与第二个子阵得到的信号矩阵进行互相关运算，得到子阵1与子阵2的互相关矩阵；将子阵1的自相关矩阵进行特征值分解得到子阵1的信号子空间矩阵，利用信号子空间矩阵的共轭转置矩阵与子阵1的自相关矩阵相乘、然后将其结果再与子阵1的信号子空间矩阵相乘，得到一个对角矩阵，将这个对角矩阵的各对角元均减去子阵1的小特征值的平均值，得到一个信号特征矩阵，再将所得信号特征矩阵中的各对角元取倒数，建立信号特征矩阵的逆矩阵；将子阵1的信号子空间矩阵的共轭转置矩阵与子阵1与子阵2的互相关矩阵相乘、其结果在再乘以子阵1的信号子空间矩阵，得到一个子阵1与子阵2互相关矩阵的过渡矩阵；该过渡矩阵与上述信号特征矩阵的逆矩阵相乘、即得到旋转不变关系矩阵；最后对旋转不变关系矩阵进行特征值分解、分解为一组特征值和与该组特征值对应的特征向量，其中所得各特征值即为不同信号来波方向的角度参数，最后经过反正弦运算，即得到不同信号的来波方向角。因而，本发明方法包括：

步骤A.划分子阵及建立各子阵的信号矩阵：

将均匀线性阵列（ULA）传感器天线的前M-1个阵元组合成子阵1，再将第2个阵元至第M个阵元组合成子阵2；将子阵1和子阵2中各阵元的接收信号分别建立接收信号矩阵、即：首先将第一时刻两子阵中各阵元得到的信号分别依次排成一列、作为对应子阵中的第一个列向量，再将下一时刻两子阵中各阵元得到的信号分别依次排列组成第二列、得到第二个列向量，以此类推按时间先后顺序将各阵元接收到的信号列向量依次排列，分别得到子阵1和子阵2的接收信号矩阵；

步骤B.建立自相关矩阵和互相关矩阵：

根据步骤A所得子阵1和子阵2两个子阵的接收信号矩阵，通过：

$R_{x_1{x}_1}=X_1{X_1}^H$

建立自相关矩阵；通过：

$R_{x_1{x}_2}=X_1{X_2}^H$

建立互相关矩阵；

上式中：

为子阵1接收信号的自相关矩阵，

为子阵1与子阵2接受信号的互相关矩阵；X₁、X₂分别为子阵1、子阵2的接收信号矩阵，(·)^H表示矩阵的共轭转置；

步骤C.建立信号子空间矩阵：

对自相关矩阵

通过：

$R_{x_1{x}_1}=\mathrm{U\Σ}{U}^H$

进行特征值分解；

其中：U为以M-1个正交的特征向量构成的酉矩阵、U=[u₁,u₂,…,u_M-1]，u₁,u₂,…,u_M-1为子阵1接收信号的自相关矩阵

的M-1个特征值对应的特征向量；∑为以M-1个特征值为对角元的对角矩阵、∑=diag(λ₁,λ₂,…λ_M-1)，λ₁,λ₂,…λ_M-1为

的M-1个特征值；(·)^H表示矩阵的共轭转置；

以不同信号的个数K为基准，将自相关矩阵

的特征值按照从大到小减值排列，利用排列后的前K个较大特征值对应的特征向量建立子阵1的信号子空间矩阵；其信号子空间矩阵为：

U_S=[u_d1,u_d2,…u_dk]

将后M-K-1个小特征值取平均值，得到λ_E待用；

上式中：u_d1，u_d2，…u_dk依次为前K个大特征值对应的特征向量；

步骤D.建立信号特征矩阵及其逆矩阵：

把子阵1接收信号的自相关矩阵

两边分别乘以步骤C所得的信号子空间矩阵U_S的共轭转置矩阵和信号子空间矩阵U_S，其结果再减去由步骤C提取待用的小特征值平均值λE构成的K×K阶的对角矩阵；得到一个K×K阶的信号特征矩阵：

B=diag(λ_d1-λ_E,λ_d2-λ_E,…λ_dK-λ_E)

其中：λ_d1,λ_d2,…λ_dK依次为自相关矩阵

特征值按照从大到小减值排列的前K个大特征值；

建立上述信号特征矩阵的逆矩阵：即取该对角矩阵中各对角元的倒数，并以各对角元的倒数建立信号特征矩阵的逆矩阵、逆矩阵为：

B^-1=diag[1/(λ_d1-λ_E),1/(λ_d2-λ_E),…,1/(λ_dK-λ_E)]

步骤E.建立过渡矩阵：

首先将由步骤C中提取待用的小特征值平均值λ_E乘以(M-1)×(M-1)阶的下三角矩阵Z得到加权矩阵：λ_EZ

其中：

利用子阵1与子阵2接收信号矩阵所得的互相关矩阵减去上述所得加权矩阵λ_EZ；其结果再两边分别乘以步骤C中所得的信号子空间矩阵U_S的共轭转置矩阵和信号子空间矩阵U_S，得到过渡矩阵

即：

$B_{x_1{x}_2}={U_S}^H(R_{x_1{x}_2}-{\mathrm{\λ}}_{E}Z)U_S$

步骤F.建立旋转不变关系矩阵：

将步骤E中得到的矩阵

左乘步骤D中得到的信号特征矩阵的逆矩阵B^-1，即得子阵1的旋转不变关系矩阵ψ，即：

$\mathrm{\ψ}=B^{-1}{B}_{x_1{x}_2};$

步骤G.获得不同信号的来波方向角：

对步骤F中所得的旋转不变关系矩阵ψ进行特征值分解，得到所需的K个特征值及与各个特征值对应的特征向量，所得特征值为：ξ₁、ξ₂…、ξ_K；将各个特征值经过自然对数和反正弦运算得到对应信号的来波方向角，即：

${\mathrm{\θ}}_k=\arcsin \left(\frac{\sqrt{-1}\mathrm{\Λ}}{2\mathrm{\πd}}\ln {\mathrm{\ξ}}_k\right)$

其中：Λ为载波波长，d为阵元间距（一般取为半波长）。

本发明采用直接建立旋转不变关系矩阵的方法，只利用矩阵的乘法和加法处理构造阵列旋转不变关系的方法有利于计算机处理、并可通过FPGA实现传感器天线阵列对单个或多个无线传播信号的来波方向进行快速测定（定向）。与现有技术比较，不但处理速度快，在相同的信噪比的情况下来波方向角的测定精度得到显提高，而且本发明方法在满足平移不变性的二维阵列中还可运用于测定空间来波方向角；因而本发明具有对设备和硬件性能要求低，对目标来波方向的测定快速、准确，精度高，流程短、运行成本低，应用范围广等特点。

附图说明

图1为本发明均匀线阵子阵划分和接收信号模型示意图；

图2为本发明具体实施方式流程图。

具体实施方式

本实施方式中无线传播信号在空间高斯白噪声的背景下、各个（来波）信号之间是互不相关的，天线阵列为均匀直线阵列，每个阵元均是全向天线，阵元个数为M=16，阵元间的间隔为d=0.083米；本实施方式中信号源为来自K=3个不同方向的远场窄带信号、且都为相同频率的载波，其波长为Λ=0.166米，且3个信号源分别以角度θ₁=20°、θ₂=30°、θ₃=40°入射到传感器均匀直线阵列；相邻两阵元中对同一信号源、其后一个阵元相比于前一个阵元的相位延迟为：φ_i=2πd sinθ_i//Λ=πsinθ_i；采样快拍数为N=1024，信噪比SNR=-15dB。

步骤A.划分子阵及建立各子阵的信号矩阵：

前15个阵元组合成子阵1，第2个到第16个阵元组合成子阵2；则两个子阵的接收数据矩阵为：

X₁=A₁S+E₁；

X₂=A₂S+E₂=A₁ΦS+E₂；

其中X₁为子阵1的接收数据矩阵，X₂为子阵2的接收数据矩阵；A₁为子阵1的方向矩阵、A₂为子阵2的方向矩阵；S为3×1024维的源信号矩阵；Φ为3×3维的旋转不变关系矩阵，E₁为子阵1的高斯白噪声矩阵，E₂为子阵2的高斯白噪声矩阵；本实施例信号与空间噪声不相关；

步骤B.建立自相关矩阵和互相关矩阵：

根据两个子阵的接收信号矩阵X₁和X₂，建立自相关矩阵

和和互相关矩阵

用阵列的导向矢量表示为如下的形式：

$R_{x_1{x}_1}=A_{1}P{A_1}^H+{\mathrm{\σ}}^{2}I$

$R_{x_1{x}_2}=A_{2}P{A_2}^H+{\mathrm{\σ}}^{2}Z$

$=A_{1}P{\mathrm{\Φ}}^H{A_1}^H+{\mathrm{\σ}}^{2}Z$

其中上标(·)^H表示矩阵的共轭转置，P=SS^H为信号自相关矩阵，σ²I和σ²Z均为15×15阶的方阵、分别为两个子阵接收噪声矩阵的自相关矩阵和互相关矩阵；

步骤C.建立信号子空间矩阵：

对自相关

进行特征值分解，得到以特征值为对角元的对角矩阵∑=diag(λ₁,λ₂,…λ₁₅)，以及咯特征向量组成的酉矩阵U=[u₁,u₂,…u₁₅]；

将自相关的特征值按照降序排列，将前3个较大特征值λ₁、λ₂、λ₃对应的特征向量u₁、u₂、u₃构成信号子空间，即信号子空间矩阵U_S=[u₁,u₂,u₃]；

将后12个较小特征值取平均值λ_E待用；

步骤D.建立信号特征矩阵及其逆矩阵：

把自相关

两边分别乘以信号子空间矩阵U_S的共轭转置矩阵和信号子空间矩阵U_S，再减去由自相关矩阵

中较小特征值的平均值λ_E构成的3×3阶的对角矩阵λ_EI；得到一个3×3阶的信号特征矩阵：

B=diag(λ₁-λ_E,λ₂-λ_E,λ₃-λ_E)

然后建立该信号特征矩阵的逆矩阵：即取该对角矩阵中各对角元的倒数，并以各对角元的倒数建立信号特征矩阵的逆矩阵、其逆矩阵为：

B^-1=diag[1/(λ₁-λ_E),1/(λ₂-λ_E),1/(λ₃-λ_E)]

步骤E.建立过渡矩阵：

将互相关

减去由步骤C中提取待用的较小特征值的均值λ_E乘以3×3阶的下三角矩阵Z得到的加权矩阵λ_EZ，得到的结果再两边分别乘以信号子空间矩阵U_S的共轭转置矩阵和信号子空间矩阵U_S，得到过渡矩阵

$B_{x_1{x}_2}={U_S}^H(R_{x_1{x}_2}-{\mathrm{\λ}}_{E}Z)U_S$

步骤F.建立旋转不变关系矩阵：

将在步骤E中得到的过渡矩阵

左乘在步骤D中得到的逆矩阵B^-1，得到的旋转不变关系矩阵ψ：

$\mathrm{\ψ}=B^{-1}{B}_{x_1{x}_2}$

步骤G.获得不同信号的来波方向角

对步骤F中所得的旋转不变关系矩阵ψ进行特征值分解，得到所需的3个特征值ξ₁、ξ₂、ξ₃及与各个特征值对应的特征向量，则3个信号源的方位角分别为：

本实施方式在信噪比SNR=-15dB时、即在信噪比较低时，其定位误差仅为时0.0167°、0.0346°及0.0746°；而在相同情况下传统的统计最优的总体最小二乘ESPRIT（TLS-ESPRIT）所测得的方位角分别为16.3254°、26.3850°、36.3465°，其定位误差依次为3.6746°、3.6150°及3.6535°；即在信噪比较低时总体最小二乘ESPRIT方法不能对信号的方位角进准确测量。

在处理速度方面，本发明随着阵元个数的增大而具有更加明显的优势；在本实施例中所采用的阵元数为16的情况下、较总体最小二乘ESPRIT方法快23.2%，当阵元数增加到100时，本发明比总体最小二乘ESPRIT方法快45.4%。

Claims (1)

1.一种对传感器天线阵列接收信号的快速处理方法，包括：

步骤A.划分子阵及建立各子阵的信号矩阵：

将均匀线性阵列传感器天线的前M-1个阵元组合成子阵1，再将第2个阵元至第M个阵元组合成子阵2；将子阵1和子阵2中各阵元的接收信号分别建立接收信号矩阵、即：首先将第一时刻两子阵中各阵元得到的信号分别依次排成一列、作为对应子阵中的第一个列向量，再将下一时刻两子阵中各阵元得到的信号分别依次排列组成第二列、得到第二个列向量，以此类推按时间先后顺序将各阵元接收到的信号列向量依次排列，分别得到子阵1和子阵2的接收信号矩阵；

步骤B.建立自相关矩阵和互相关矩阵：

根据步骤A所得子阵1和子阵2两个子阵的接收信号矩阵，通过：

$R_{x_1{x}_1}=X_1{X_1}^H$

建立自相关矩阵；通过：

$R_{x_1{x}_2}=X_1{X_2}^H$

建立互相关矩阵；

上式中：

为子阵1接收信号的自相关矩阵，

为子阵1与子阵2接收信号的互相关矩阵；X₁、X₂分别为子阵1、子阵2的接收信号矩阵，(·)^H表示矩阵的共轭转置；

步骤C.建立信号子空间矩阵：

对自相关矩阵通过：

$R_{x_1{x}_1}=\mathrm{U\Σ}{U}^H$

进行特征值分解；

其中：U为以M-1个正交的特征向量构成的酉矩阵：U=[u₁,u₂,…,u_M-1]，u₁,u₂,…,u_M-1为子阵1接收信号的自相关矩阵的M-1个特征值对应的特征向量；Σ为以M-1个特征值为对角元的对角矩阵、Σ=diag(λ₁,λ₂,…λ_M-1)，λ₁,λ₂,…λ_M-1为

的M-1个特征值；(·)^H表示矩阵的共轭转置；

以不同信号的个数K为基准，将自相关矩阵

的特征值按照从大到小减值排列，利用排列后的前K个较大特征值对应的特征向量建立子阵1的信号子空间矩阵；其信号子空间矩阵为：

U_S=[u_d1,u_d2,…u_dk]

将后M-K-1个小特征值取平均值，得到λ_E待用；

上式中：u_d1,u_d2,…u_dk依次为前K个大特征值对应的特征向量；

步骤D.建立信号特征矩阵及其逆矩阵：

把子阵1接收信号的自相关矩阵

两边分别乘以步骤C所得的信号子空间矩阵U_S的共轭转置矩阵和信号子空间矩阵U_S，其结果再减去由步骤C提取待用的小特征值平均值λ_E构成的K×K阶的对角矩阵；得到一个K×K阶的信号特征矩阵：

B=diag(λ_d1-λ_E,λ_d2-λ_E,…λ_dK-λ_E)

其中：λ_d1,λ_d2,…λ_dK依次为自相关矩阵特征值按照从大到小减值排列的前K个大特征值；

建立上述信号特征矩阵的逆矩阵：即取该对角矩阵中各对角元的倒数，并以各对角元的倒数建立信号特征矩阵的逆矩阵，该逆矩阵为：

B^-1=diag[1/(λ_d1-λ_E),1/(λ_d2-λ_E),…,1/(λ_dK-λ_E)]

步骤E.建立过渡矩阵：

首先将由步骤C中提取待用的小特征值平均值λ_E乘以(M-1)×(M-1)阶的下三角矩阵Z得到加权矩阵：λ_EZ

其中：

利用子阵1与子阵2接收信号的互相关矩阵

减去上述所得加权矩阵λ_EZ；其结果再两边分别乘以步骤C中所得的信号子空间矩阵U_S的共轭转置矩阵和信号子空间矩阵U_S，得到过渡矩阵

即：

$B_{x_1{x}_2}={U_s}^H(R_{x_1{x}_2}-{\mathrm{\λ}}_{E}Z)U_s$

步骤F.建立旋转不变关系矩阵：

将步骤E中得到的矩阵左乘步骤D中得到的信号特征矩阵的逆矩阵B^-1，即得子阵1的旋转不变关系矩阵ψ，即：

$\mathrm{\ψ}=B^{-1}{B}_{x_1{x}_2};$

步骤G.获得不同信号的来波方向角：

对步骤F中所得的旋转不变关系矩阵ψ进行特征值分解，得到所需的K个特征值及与各个特征值对应的特征向量，所得特征值为：ξ₁、ξ₂…、ξ_K；将各个特征值经过自然对数和反正弦运算得到对应信号的来波方向角，即：

${\mathrm{\θ}}_k=\arcsin \left(\frac{\sqrt{-1}\mathrm{\Λ}}{2\mathrm{\πd}}1n{\mathrm{\ξ}}_k\right)$

其中：Λ为载波波长，d为阵元间距。

Images