互质阵列中基于迭代稀疏重构的DOA估计方法

摘要

本发明公开了互质阵列中一种基于迭代稀疏重构的DOA估计方法，接收天线阵列采用非线性互质阵列，通过对接收信号的二阶统计特性协‑方差矩阵的向量化处理，可确定孔径长度更大的差分阵列，提高检测能力。对目标所在的角度域进行离散化处理，目标可视为稀疏分布于栅格点之上或附近，并构建出对数和形式的稀疏信号重建问题。利用对数和函数的凸紧上界，重建该原始稀疏问题，以迭代的方式动态调整更新角度域的离散点以逼近目标真实的到达角度。

说明书

技术领域

本发明涉及通信信号处理技术领域，尤其涉及一种使用互质阵列的基于迭代稀疏重构的DOA估计方法。

背景技术

基于天线阵列的信号处理方法被广泛应用于无线通信、电磁场、雷达、声呐等诸多领域。波达方向(DOA)估计是阵列信号处理领域中的一个重要问题。

传统的估计方法通常研究阵元间隔为半波长的均匀线性阵列，且适用于检测目标数少于阵元个数的场合，例如对N天线均匀线性阵列而言，传统估计方法(如基于子空间的估计方法等)至多可检测N-1个目标，非线性互质阵列通过利用协方差矩阵的特性，将原始阵列构建成具有更多虚拟天线和更大孔径长度的差分阵列，进而可显著提高其自由度，即检测能力，传统估计方法通常需要目标个数的先验信息和大样本才能估计出目标的波达方向，在小样本及目标个数未知条件下不适用，此外，传统估计方法也难以适用于信噪比低的情况，即可能无法检测到发射功率较低(例如小于噪声功率)的弱目标。

近年来，随着稀疏重构理论的发展，稀疏表示已逐渐开始应用于小波去噪、雷达成像等领域，基于稀疏信号重构的DOA估计方法能够充分发掘互质阵列自由度高的优势，可检测目标数目较传统方法显著提高，然而，常规的稀疏重构方法需要将角度域进行静态的栅格化处理，建立以l₀范数(通常以l_p范数近似，p≤1)最小化为目标的稀疏优化问题，这类方法的问题在于估计精度严重依赖于初始的栅格化过程，若目标大部分位于或极其靠近栅格点，则估计性能良好；反之，若大部分目标流离于栅格点之外，则估计性能难以保证，因此，由于角度域的静态栅格化处理造成的栅格失配问题会严重影响重构的效果。

发明内容

本发明针对上述现有技术存在的不足，提出了一种适用于互质阵列的基于迭代稀疏重构的DOA估计方法，本发明的方法具有更高的自由度，可检测比实际阵元数更多的目标，本发明的方法的估计性能不依赖于初始的栅格化处理，具有更强的灵活性，适用于小样本、低信噪比及目标个数未知等复杂环境，其中所述DOA全称为Direction of Arrival,即信号的波达方向。

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

一种互质阵列中基于迭代稀疏重构的DOA估计方法，它包括以下步骤：

(1)建立接收阵列模型

用非线性的互质阵列作为接收阵列采集信号，其中互质阵列可分解为两个均匀线性子阵列，每个子阵列的阵元间隔均大于半波长；

(2)构建虚拟差分阵列

对互质阵列接收信号的协方差矩阵进行向量化处理，在不同阵元的差分位置形成虚拟阵元，可将原始阵列等效成具有更多虚拟阵元及更大孔径尺度的差分阵列；

(3)栅格化处理，建立稀疏优化问题

将目标区域进行栅格化处理，分析差分阵列中的接收信号，建立以对数和函数为优化目标的稀疏优化问题，兼顾估计结果的稀疏性和准确性；

(4)DOA估计迭代实现

利用对数和函数的凸紧上界重构原始稀疏问题，以迭代方式逐步更新栅格点位置来使目标函数最小化，逐渐逼近真实信号源的方向，直至满足终止条件，所述方法可动态调整栅格点的位置，能有效克服角度域的静态离散化造成的栅格失配问题。

本发明的有益效果为：

1、本发明利用互质阵列构建出具有更多阵元数目和更大孔径尺寸的虚拟差分阵列，显著提高了自由度，可实现多于阵元数目的目标检测；

2、相对于传统角度域的静态栅格化，本发明以迭代方式动态调整栅格点来逐渐逼近目标的真实位置，可避免估计性能对栅格初始位置的依赖，改善由于静态的栅格化造成的失配问题，提高估计精度和分辨率；

3、本发明无需目标个数的先验信息及对协方差矩阵的满秩性要求，可适用于目标个数未知及小样本环境，具有很强的灵活性。在单样本条件下，该方法可实现目标的动态追踪；

4、本发明鲁棒性好，对信噪比要求低，可实现低发射功率(小于噪声功率)目标的检测。

附图说明

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

图2为本发明互质阵列及其子阵列示意图；

图3为本发明的仿真实验中归一化功率谱图；

图4为本发明的仿真实验中估计精度图a；

图5为本发明的仿真实验中估计精度图b。

具体实施方式

下面结合附图对本发明作进一步的详细描述：

本发明方法的流程图如图1所示，具体实现过程如下：

(1)构建非线性互质阵列模型，获取接收信号；

(2)计算接收信号的协方差矩阵并进行向量化处理，建立虚拟差分阵列；

(3)对角度域进行栅格化处理，建立稀疏优化问题；

(4)以迭代方式动态调整栅格点位置，直至满足终止条件。分析稀疏解，确定最终波达方向。

本发明的DOA估计方法过程如下：

1.互质阵列及其接收信号

本发明所涉及的互质阵列如图2所示，该阵列可分解为两个均匀线性子阵列，其中子阵列1包含M₁个天线，相邻天线间隔为M₂λ/2，子阵列2包含2M₂个天线，相邻天线间隔为M₁λ/2，这里M₁和M₂为互质的正整数，λ表示载波的波长，子阵列1和2的全体构成非线性互质阵列，由于子阵列1和2共享第一个天线，因此互质阵列的天线数为M＝M₁+2M₂-1。

有未知个数(假定为K)的不相关目标从不同方向Θ＝[θ₁,θ₂,…,θ_K]到达互质阵列，该阵列在时刻t(1≤t≤T)的接收信号为

$x\left(t\right)=\underset{k=1}{\overset{K}{\Sigma }}a\left({\theta }_k\right)s_k\left(t\right)+n\left(t\right)=As\left(t\right)+n\left(t\right)---\left(1\right)$

其中，A＝[a(θ₁),a(θ₂),…,a(θ_K)]表示由阵列位置确定的已知的阵列流形矩阵，s(t)＝[s₁(t),s₂(t),…,s_K(t)]^T表示K个目标的发射信号矢量，n(t)为独立同分布加性高斯白噪声矢量，上标T表示转置。

2.协方差矩阵及虚拟差分阵列

接收信号x(t)的协方差矩阵可表征为

$R_{xx}=E\left[x\left(t\right)x^H\left(t\right)\right]=\underset{k=1}{\overset{K}{\Sigma }}{\sigma }_k^{2}a\left({\theta }_k\right)a^H\left({\theta }_k\right)+{\sigma }^2{I}_{M_1+2M_2-1}---\left(2\right)$

其中，和σ²分别表示第k个信号的功率和噪声功率，E表示期望，上标H表示共轭转置，>1+2M₂-1维单位阵，矩阵R_xx的第(m,n)项为可视为在l_m-l_n处存在的虚拟天线的接收信号，在R_xx中，由第m个和第n个天线的差值产生的虚拟天线在位置l_m-l_n(1≤m,n≤M₁+2M₂-1)，l_m和l_n分别表示第m和第n个天线的实际位置。

对R_xx进行向量化处理，则有

z＝vec(R_xx)＝Φ(θ₁,θ₂,…,θ_K)p+σ²1_n,>

其中和vec表示向量化处理，表示Kronecker乘积，z为虚拟差分阵列的接收信号，Φ(θ₁,θ₂,…,θ_K)表示虚拟差分阵列的阵列流行矩阵。

3.栅格化处理，建立稀疏优化问题

为利用稀疏方法进行DOA估计，需要将目标角度域进行栅格化， 因此可建立稀疏优化问题为

$\begin{array}{ccc}\underset{{\theta }^g,p,{\sigma }^2}{min}\left|\right|p|{|}_0,& s.t.& z=\Phi \left({\theta }^g\right)p+{\sigma }^2{1}_n\end{array}---\left(4\right)$

其中||·||₀表示0-范数，该优化问题表示的含义为在给定静态栅格点，即在给定虚拟阵列流行矩阵的前提下，用尽可能少的信号功率p去重构接收信号z，该问题是NP-难问题，计算量极大，

为此，利用对数和函数去近似(4)中的0-范数，并构建无约束优化问题为

$\begin{array}{cc}\underset{{\theta }^g,p,{\sigma }^2}{min}& \eta \underset{i=1}{\overset{D}{\Sigma }}log(p_i^2+ϵ)+\frac{1}{2}\left|\right|z-\Phi \left({\theta }^g\right)p-{\sigma }^2{1}_n|{|}_2^2\end{array}---\left(5\right)$

其中ε>0用来确定对数函数的存在性，||·||₂表示最小二乘代价函数，η>0衡量稀疏性和最小二乘代价之间的均衡，由于对数函数的非凸性，优化问题(5)容易陷入局部最优解，进一步，利用对数和函数的凸上界函数

$Q\left(p\right|{\hat{p}}^{\left(t\right)})=\underset{i=1}{\overset{N}{\Sigma }}(\frac{|p_i{|}^2+ϵ}{|{\hat{p}}_i^{\left(t\right)}{|}^2+ϵ}+log(|{\hat{p}}_i^{\left(t\right)}{|}^2+ϵ)-1),---\left(6\right)$

来代替优化问题中的对数和函数，式中，p_i为p中的第i个元素，为p_i在第t次迭代中的估计，替代中忽略掉与未知变量{p,θ^g,σ²}无关的项之后，优化问题(5)可转化为

$\begin{array}{cc}\underset{{\theta }^g,p,{\sigma }^2}{min}& {\mathrm{\eta p}}^H{D}^{\left(t\right)}p+\frac{1}{2}\left|\right|z-\Phi \left({\theta }^g\right)p-{\sigma }^2{1}_n|{|}_2^2,\end{array}---\left(7\right)$

其中

4.DOA估计迭代实现

该迭代方法的具体实施步骤如下：

步骤1：初始化离散角度集合θ^g,(0)，对应的信号功率集合和噪声功率σ²_,⁽⁰⁾,并令t＝1，

步骤2：根据当前角度集合θ^g,(t-1)和噪声功率σ^2,(t-1)，优化问题(7)对p求导并置零，计算得

${\hat{p}}^{\left(t\right)}={({\Phi }^H\left({\theta }^{g,(t-1)}\right)\Phi \left({\theta }^{g,(t-1)}\right)+2{\mathrm{\eta D}}^{(t-1)})}^{-1}{\Phi }^H\left({\theta }^{g,(t-1)}\right)(z-{\sigma }^{2,(t-1)}{1}_n)---\left(8\right)$

步骤3：根据当前角度集合θ^g_,^(t-1)和信号功率估计计算噪声功率σ²_,^(t)为

${\sigma }^{2,\left(t\right)}=\frac{1}{2N}(z^H{1}_n+1_n^{H}z-{\left({\Phi }^{(t-1)}{\hat{p}}^{\left(t\right)}\right)}^H{1}_n-1_n^H\left({\Phi }^{(t-1)}{\hat{p}}^{\left(t\right)}\right))---\left(9\right)$

步骤4：根据当前信号功率估计构建对数和函数的凸上界函数更新

步骤5：将估计值σ²_,^(t)和D^(t)代入优化问题(7)，优化问题变为>

$\begin{array}{cc}\underset{{\theta }^g}{min}& f\left({\theta }^g\right)=-{(z-{\sigma }^2{1}_n)}^H\Phi \left({\theta }^g\right){({\Phi }^H\left({\theta }^g\right)\Phi \left({\theta }^g\right)+2{\mathrm{\eta D}}^{\left(t\right)})}^{-1}{\Phi }^H\left({\theta }^g\right)(z-{\sigma }^{2,\left(t\right)}{1}_n)\end{array}---\left(10\right)$

由于Φ(θ^g)是关于θ^g的非线性函数，直接获取最优θ^g难以实现，可采用迭代方式逐步靠近最优的θ^g，为此寻找新的估计θ^g,(t)满足下式

f(θ^g,(t))≤f(θ^g,(t-1))

θ^g_,^(t)可根据梯度下降法估计为

${\theta }^{g,\left(t\right)}={\theta }^g-\mu \frac{\partial f\left({\theta }^g\right)}{\partial {\theta }^g}{|}_{{\theta }^g={\theta }^{g,(t-1)}}$

其中μ为较小的正数，

令t＝t+1，

步骤6：若满足终止条件，算法结束，否则跳至步骤2。

以下通过仿真对比本发明方法和其他传统方法(如空间平滑方法)说明本发明的优越性能：

本发明仿真实验采用M₁＝5和M₂＝3的互质阵列模型，其天线总个数为10，存在K＝11个等功率信号源，其角度分别为[-49.3，-37.2，-26.8，-17.3，-8.3，0.45，9.2，18.3，27.8，38.3，50.6]度，将角度域按照3度的间隔进行等间隔栅格化，信噪比定义为输入功率与噪声功率之比。

图3为本发明的仿真实验中归一化功率谱图，其中图2中的横坐标表示波达方向，纵坐标表示归一化能量，虚线表示真实角度，上图实线表示空间平滑算法估计的角度，下图实线表示本发明方法估计的角度，如图3所示，本发明方法能成功检测所有的目标，而空间平滑方法漏掉了其中一个，因此本发明方法具有更强的检测能力。

图4为本发明的仿真实验中估计精度图，定量分析了本发明的估计精度，图4中横坐标表示信噪比，纵坐标表示估计均方误差，图5中横坐标表示样本数目，纵坐标表示估计的均方误差，由图4可见，在不同信噪比条件下，本发明的估计均方误差均小于空间平滑算法的估计误差，在低信噪比环境下更加明显，由图5可见，在不同样本数目条件下，本发明的估计均方误差也均小于空间平滑算法的估计误差，在小样本条件下更加明显，显然，本发明的DOA估计精度高于现有的空间平滑算法。