一种利用互质阵的波达方向估计方法与流程

本发明涉及阵列信号处理领域，尤其是一种波达方向估计方法。

背景技术：

波达方向估计在雷达、声纳和无线通信中，有很广泛的应用，然而传统的波达方向估计方法只能求解目标数少于阵元数的情况，因此如何用少量的阵元检测更多的目标是一个值得研究的问题。实际中，一般常用的N元均匀线列阵，最多只能估计出N-1个目标的来波方向。近几年，提出了一种新型几何结构的线列阵——互质阵，使得估计的目标数目远超过阵元数目N，事实上，一个有N+M阵元的互质阵所能估计的目标数目可以达到O(MN)，即与MN同阶大小的目标数目，由于互质阵阵元位置的巧妙分布，在经过数学运算处理后，可以形成孔径更大的虚拟阵。具体来说，上述提到的N+2M-1个物理阵元的互质阵可以虚拟出孔径为2MN+1的均匀线列阵。

对于互质阵，已经提出了一种基于互质阵的空域平滑MUSIC算法，对于使用N+2M-1个物理阵元的互质阵，使用该方法最多可以估计出MN个目标，尽管这个算法估计出的目标数目已经远大于同等数目物理阵元的均匀线列阵，但是由于使用了空域平滑方法，因此导致虚拟阵元减少到MN+1，检测性能降低。最近几年，由D.Donoho、E.Candes及华裔科学家T.Tao等人提出了一种新的信息获取指导理论，即压缩感知，该理论一经提出，就在信息论、信号处理、模式识别、无线通信等领域受到高度关注，该理论认为，对于一个稀疏信号或者是可以稀疏表示的信号，可以对该信号进行少量观测后，然后采用特定的重构算法重构原始信号，对于信号x的有限观测y，满足y＝Ax，式中观测矩阵M＜N，若满足一定条件，便可从有限观测y中重构原始的稀疏信号x。

技术实现要素：

为了克服现有技术的不足，本发明基于压缩感知中的稀疏重构思想，对于互质阵构造基于OMP算法的波达方向估计方法，克服了空域平滑MUSIC算法导致虚拟阵元减少的问题，使得检测目标数目超过MN，并且仿真证明该方法在高信噪比时有很好的估计精度。

本发明解决其技术问题所采用的技术方案具体包括如下步骤：

第一步：互质阵结构设计

互质阵由两个间距分别为Md和Nd的均匀线列阵嵌套组成，其中M和N为两个互为质数的常数，d表示接收信号的半波长，一个均匀线列阵有N个阵元，阵元间距为Md，另外一个均匀线列阵有2M个阵元，阵元间距为Nd，两个均匀线列阵的第一个阵元重合，共有N+2M-1个阵元，由于M和N都是质数，所以称这种线列阵为互质阵，互质阵是一种特殊结构的非均匀线列阵；

第二步：估计互质阵接收信号的协方差矩阵

互质阵接收信号的协方差矩阵R_xx通过公式(1)估计：

式(1)中，表示互质阵接收信号向量，T表示阵列快拍数；

第三步：构造虚拟均匀线列阵的接收向量y

构造大小为(N+2M-1)×(N+2M-1)标识矩阵B，B中第i行第j列的元素b_ij＝i-j，标识矩阵B和协方差矩阵R_xx分别按列拉直，得到标识向量从标识向量b中按照从小到大的顺序依次寻找元素值为-MN到MN的位置，并按顺序记录下2MN+1个位置信息，然后从向量z中依次将对应位置的元素提取出来，构造出虚拟均匀线列阵的接收向量

第四步：待检测的角度区域剖分网格，构造阵列流形矩阵A

将待检测的角度区域离散化剖分网格，形成网格向量θ＝[θ₁,θ₂,…,θ_D]^T，其中θ_k表示离散的网格角度，k＝1,…,D，D表示网格数目，且使得剖分的网格数目D大于信号数目，由公式(2)和公式(3)构造阵列流形矩阵A：

A(θ₁,θ₂,…,θ_D)＝[a(θ₁),a(θ₂),…a(θ_k),...,a(θ_D)] (2)

式(2)中，a(θ_k)为对应阵列的阵列流形向量，k＝1，2，...，D，式(3)中，j表示虚数单位，d表示接收信号的半波长，λ表示接收信号的波长，θ_k表示离散的网格角度；

第五步：利用OMP算法估计稀疏信号

阵列流形矩阵A(θ₁,θ₂,…,θ_D)和虚拟均匀线列阵的接收向量y满足如下方程：

式(4)中p为对应于剖分网格θ＝[θ₁,θ₂,…,θ_D]^T的稀疏向量，表示噪声向量，稀疏向量p采用压缩感知中的OMP算法求解可得稀疏向量估计值

第六步：由稀疏向量得到信号波达方向

若稀疏向量估计值中的第i项非0，则表示对应的θ_i方向具有信号，否则表示没有信号。

本发明相比传统的波达方向估计采用均匀线列阵，由于采用了互质阵这种新型的嵌套阵，使得波达方向估计性能显著的提高了估计的目标数量，在阵列多目标估计时具有很好的估计性能；在同样使用互质阵时，该方法相比于空间平滑的MUSIC算法而言，由于使用了压缩感知的OMP算法，使得估计的目标数目超过NM个目标，在同样阵元数量下，进一步提高了互质阵的估计目标数目；由于使用了OMP算法，可以简单高效的对稀疏信号进行重构，最终快速高效的给出波达方向估计结果；在高信噪比时，本发明方法的估计精度优于空间平滑的MUSIC算法，总体而言，本发明的性能不低于空间平滑的MUSIC算法。

附图说明

图1是本发明进行波达方向估计的方法流程图。

图2是本发明互质阵的几何结构。

图3是本发明与空间平滑MUSIC算法的波达方向估计结果。

图4是本发明对16个来波信号的波达方向估计结果。

图5是本发明与空间平滑MUSIC算法的单目标的估计精度对比结果。

图5中MUSIC表示空间平滑MUSIC算法，OMP表示本发明提出的基于OMP算法的波达方向估计方法。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

第一步：互质阵结构设计

选取两个互为质数的常数M和N，用常数d表示接收信号的半波长，互质阵由两个间距分别为Md和Nd的均匀线列阵嵌套组成，其中一个均匀线阵有N个阵元，阵元间距为Md，另外一个线阵有2M个阵元，阵元间距为Nd，两个线列阵第一个阵元重合，总共有N+2M-1个阵元，由于M和N都是质数，所以称这种线列阵为互质阵，互质阵的结构如图1所示，互质阵是一种特殊结构的非均匀线列阵，因此可以在线列阵的模型上进一步推导；

首先考虑包含N个阵元的线阵，阵元位置为l_i，假设K个目标方位为并且设K个目标信号为窄带信号，那么阵列接收信号表示为：

式中1≤t≤T，s_k(t)是第k个来波信号，n(t)是独立同分布的白噪声，是阵列流形向量，且

的第i个元素代表第k个信号在第i个阵元的时间延迟所带来的相位变化。

第二步：估计互质阵接收信号的协方差矩阵

互质阵接收信号的协方差矩阵R_xx通过公式(1)估计：

式中，表示互质阵接收信号向量，T表示阵列快拍数；

设信号s_k(t)服从方差为的独立高斯分布，考虑每个阵元接收数据的二阶统计量，求解阵列接收信号x(t)的协方差矩阵R_xx

式中σ²是噪声功率，将矩阵R_xx按列拉直，根据矩阵Kronecker积(也称直积)的运算性质，由式(7)得到向量z：

其中，表示第k个发射信号的方差，表达式如下：

代表除了第i个位置为1其余元素都为0的列向量，根据(8)式将向量z看作是信源向量q经过矩阵Φ观测后接收到的信号，矩阵中存在一个阵列流形矩阵，并且这个矩阵是一个具有更多虚拟阵元的线列阵的阵列流形矩阵。

第三步：构造虚拟均匀线列阵的接收向量y

构造大小为(N+2M-1)×(N+2M-1)标识矩阵B，B中第i行第j列的元素b_ij＝i-j，标识矩阵B和协方差矩阵R_xx分别按列拉直，得到标识向量从标识向量b中从小到大的顺序依次寻找元素值为-MN到MN的位置，按顺序记录下共2MN+1个位置信息，假设b＝[-3,-1,1,5,0]，从b中以小到大的顺序依次寻找元素值为-1到1的位置，记录下的对应位置信息就是[2,5,3]，然后从向量z中依次将对应位置的元素提取出来，构造出虚拟均匀线列阵的接收向量

第四步：待检测的角度区域剖分网格，构造阵列流形矩阵A

将待检测的角度区域离散化剖分网格，形成网格向量θ＝[θ₁,θ₂,…,θ_D]^T，其中θ_k(k＝1,…,D)表示离散的网格角度，D表示网格数目，且使得剖分的网格数目D大于信号数目，取D为信号数目的10倍以上，由下面公式构造阵列流形矩阵A：

A(θ₁,θ₂,…,θ_D)＝[a(θ₁),a(θ₂),…,a(θ_D)] (2)

式(2)中，a(θ_k)为对应阵列的阵列流形向量，k＝1，2，...，D，式(3)中，j表示虚数单位，d表示接收信号的半波长，λ表示接收信号的波长，θ_k表示离散的网格角度；

第五步：利用OMP算法估计稀疏信号

阵列流形矩阵A(θ₁,θ₂,…,θ_D)和虚拟均匀线列阵的接收向量y满足如下方程：

式(4)中p为对应于剖分网格θ＝[θ₁,θ₂,…,θ_D]^T的稀疏向量，表示噪声向量，稀疏向量p采用压缩感知中的OMP算法求解可得稀疏向量估计值

第六步：由稀疏向量得到信号波达方向

若稀疏向量估计值中的第i项非0，则表示θ_i方向具有信号，否则表示没有信号。

在本发明中，更具体的假设这个线列阵为如图1所示的互质阵，由两个间距分别为Md和Nd的均匀线列阵嵌套组成，一共有N+2M-1个物理阵元，其中一个均匀线阵有N个阵元，另外一个线阵有2M个阵元，两个线阵的第一个阵元重合。假设信号的协方差矩阵R_xx已知，矩阵大小为(N+2M-1)×(N+2M-1)，根据公式(4)，由直积的性质可以知道，矩阵Φ共有(N+2M-1)²个行向量，已经证明了，如果N,M都是质数，这些行向量的子集对应于一个孔径更大的均匀线列阵的阵列流形矩阵A，其阵元位置从-MNd到MNd，阵元间距为d，这时虚拟均匀线列阵的孔径可以达到(2MN+1)d，下面由矩阵构造均匀线列阵的阵列流形矩阵

从Φ中选择相应的行向量，并对这些行向量进行排序得到A，使得阵列流形矩阵A中的第(n,k)个元素为k＝1,2,…,K,n＝-MN,…,0,…,MN，最终形成大小为(2NM+1)×K的矩阵A，此时A便对于一个孔径为(2MN+1)d的均匀线列阵的阵列流形矩阵，对于式(8)左端的向量z，按照步骤3中所述构造向量y的方法选择相应的项并且排序得到向量y，y便是对应于此虚拟均匀线列阵的接收向量，那么A、y满足

式中向量是向量1_n去掉相应项排序得到的向量，并且第NM+1项为1，其余均为0。

下面对压缩感知的稀疏重构思想进行波达方向估计进行说明：

首先对待检测的角度区域离散化：θ＝[θ₁,θ₂,…,θ_D]^T，且D＞＞K，即剖分的网格数目大于信号数目，构造阵列流形矩阵A(θ₁,θ₂,…,θ_D)：

A(θ₁,θ₂,…,θ_D)＝[a(θ₁),a(θ₂),…,a(θ_D)] (2)

式中是一个阵元数为2NM+1的均匀线列阵的阵列流形向量。因此可以将公式(4)改写为

其中，y是根据互质阵的接收信号x(t)计算得到的虚拟均匀线列阵的接收向量，如果第k个信号的来波方向为θ_i，那么向量p的第i个位置假设所有目标的位置恰好位于剖分的网格点上，那么只在K个来波方向上是非零的，此时，p为K稀疏向量。对于公式(9)，如果给定y和A，通过压缩感知中的稀疏重构算法来求解方程。

根据压缩感知理论，需要求得满足公式(10)的最稀疏解因此稀疏重构问题表示为：

式中||·||₀表示向量非零项的数目，为了表示方便A＝A(θ₁,θ₂,…,θ_D)，假设噪声方差σ²是未知的，为了求解式(10)，本发明采用OMP算法求解，OMP算法的迭代步骤如下：

输入：阵列流形矩阵A，虚拟均匀线列阵的接收向量y，目标信号数量K；

输出：p的K稀疏逼近误差向量E；

初始化：余量E⁰＝y，重构信号p⁰＝0，索引集迭代次数n＝0；

步骤1：计算余量和构造阵列流形矩阵A(θ₁,θ₂,…,θ_D)的每一列之间的内积Gⁿ＝A^TE^n-1，E^n-1是第n次迭代的余量，n＝1时，E^n-1＝E⁰；

步骤2：找出Gⁿ中绝对值最大的元素对应的位置，即

步骤3：更新索引集Γⁿ＝Γ^n-1∪{k}及原子集合

步骤4：利用最小二乘法求得近似解

步骤5：更新余量

步骤6：判断迭代是否满足停止条件，即迭代次数n是否大于K，当迭代次数大于K则满足停止条件，若满足停止条件，则E＝Eⁿ，输出否则令n←n+1，转步骤1，其中，是的第i项。

以下对本发明的方法进一步描述：在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程。

考虑包含10个物理阵元的互质阵，即对图2中的互质阵取N＝5，M＝3，具体来说，第一层阵元位置在[0,3,6,9,12]d，第二层阵元位置在[0,5,10,15,20,25]d，取d为半波长值，选取若干窄带信号，来波方向均匀分布在-60度到60度的区间内，频率均取为f＝1000Hz，采样频率fs＝8192Hz，快拍数500，信噪比为0dB。

图3中选取了15个不同方向的窄带信号仿真两种算法的检测性能。图4中左图是按照空间平滑矩阵MUSIC算法得到的结果。图3中右图展示了基于OMP算法的波达方向估计。首先从-60度到60度以0.5度为步长划分网格，即D＝240，然后根据式(1)到(3)，计算构造式(4)的等式方程，最后执行OMP算法，得到波达方向估计结果。图4中虚线表示原始信号方向，实线表示估计结果。可以看出，两种算法均能得到很好的估计效果，然而空间平滑MUSIC算法最多只能估计出MN＝15个信号的来波方向。对于上述提到的10元互质阵，当有16个来波信号时，空间平滑MUSIC算法便会出现问题，MUSIC算法最多只能估计出NM＝15个信号，当16个信号时，方法无解，若依然按照15个信号的输入求解，估计结果与真实值相差很大，此时基于OMP算法的波达方向估计将展现出优势。

图4给出了16个来波信号时，使用OMP算法进行波达方向估计的结果，可以看出OMP算法的波达方向估计结果能给出准确的估计结果，而这种结果是空间平滑MUSIC算法不能达到的。

在图5给出的仿真图中，本发明考虑单个目标估计时，两种算法的估计精度随着信噪比的变化情况。不同于前面的仿真，这次在-90度到90度的范围内，随机产生一个信号方向，横轴表示的信噪比范围以2dB为步长，从-20dB到10dB；纵轴表示的误差采用估计角度与真实角度差值的绝对值。每个点执行5000次蒙特卡洛实验，其余条件与前面仿真相同。从图5中的结果可以看出，对于单目标的波达方向估计，在低信噪比时，两种算法的估计性能基本一样，但是在高信噪比时，OMP算法的估计精度优于空间平滑的MUSIC算法。