CN104459667B - 一种基于clean的稀疏阵列波达方向doa估计方法 - Google Patents

Abstract

本发明涉及一种基于CLEAN的稀疏阵列波达方向DOA估计方法，包括以下步骤：利用最小描述长度(MDL)准则来计算目标个数；计算所有快拍数据形成波束的幅度平方累加结果；在计算结果中查找最大值对应的波束号；利用等信号测角法计算“精确波束号”；计算目标角度；判断所有信号是否均已检测；应用CLEAN思想清除检测到的目标响应，转至进行迭代。该方法能够很好的抑制强目标与高旁瓣作用结果对弱目标角度估计的影响，达到很高的测角精度，同时，该方法耗时很短且耗时对阵元数不敏感。

Description

一种基于CLEAN的稀疏阵列波达方向DOA估计方法

技术领域

本发明属于雷达信号处理技术领域，具体涉及一种基于CLEAN的稀疏阵列波达方向DOA估计方法。

背景技术

在雷达和通信***中，阵列天线应用极其广泛。到目前为止，均匀间隔分布的天线是最简单且应用最广的阵列天线之一。然而均匀间隔分布的天线阵列存在出现栅瓣的可能性，为避免栅瓣的出现，通常要求天线阵列的间距不大于波长的一半。因此，如果想要获得较高的角度分辨力就需要很多的阵元，这不仅会增加天线的成本，而且产生的大量数据也会增加数字信号处理***的负担。

对均匀分布的阵列进行稀疏化处理，也就是从均匀阵列中随机去掉一些阵元，使得阵元不再规则排列得到稀疏阵列。稀疏阵列能够有效的抑制栅瓣的出现，并且能够提高角度分辨力，但它的峰值旁瓣和均值旁瓣都比较高，强目标与高旁瓣的相互作用的响应有可能会大于弱目标的响应，有出现假目标的可能性。

波达方向(DOA，Direction of Arrival)估计一直是阵列信号处理领域的热点之一，精确的DOA估计对后续处理至关重要。当前对DOA估计算法的研究主要集中在MUSIC和ESPRIT两种出色的算法之上。

MUSIC(Multiple Signal Classification)算法具有很好的分辨力，但该算法需要在所有可能的角度进行谱峰搜索，运算量极大，而且测角精度和所选的搜索角度的间隔相关。Root-MUSIC是在MUSIC的基础上提出的，其无需进行谱峰搜索，直接得到闭式解。不过，Root-MUSIC算法的策略是求方程的根，随着阵元数的增加，方程求根的时间会急剧增加，所以在阵元数较多的场合Root-MUSIC会增加时间上的负担。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种基于CLEAN的稀疏阵列波达方向DOA估计方法，解决强弱目标共存时，强目标与稀疏阵列高旁瓣相作用产生的假目标对应的响应大于弱目标响应的问题。

技术方案

一种基于CLEAN的稀疏阵列波达方向DOA估计方法，其特征在于步骤如下：

步骤1：利用最小描述长度MDL准则计算CLEAN的稀疏阵列信号中的目标个数；

步骤2：T>0时，计算所有快拍数据形成波束的幅度平方累加结果否则结束计算过程；其中：y_s为第s个快拍数据x_s利用FFT形成的波束数据，y_s＝FFT(x_s)；L为快拍数；

步骤3：在幅度平方累加结果P中查找最大值对应的波束号；

步骤4：采用等信号测角法计算“精确波束号”p_t，

$p_t=\frac{(p_1+p_2)}{2}+k\·\mathrm{sign}\·\min (\frac{|P\left(p_1\right)-P\left(p_2\right)|}{P\left(p_1\right)+P\left(p_2\right)|},0.5)$

其中，p₁为P中最大值对应的波束号，p₂为次大值的波束号，k为阵列对应的比幅测角斜率，sign为符号位；

$p_2=\left\{\begin{array}{cc}{p}_1-1,& P(p_1-1)>P(p_1+1)\\ p_1+1,& P(p_1-1)\≤P(p_1+1)\end{array}\right.$

所述 $\mathrm{sign}=\left\{\begin{array}{cc}-1,& p_1<p_2\\ 1,& p_1\&GreaterEqual;p_2\end{array}\right.$

步骤5：计算目标角度，

其中，λ为波长，d为阵元间距，M为FFT点数；

步骤6：应用CLEAN思想清除检测到的目标响应，同时将T值减一，转至步骤2进行迭代。

所述稀疏阵列是在均匀间隔阵列中随机去除一些阵元形成的，用数组A_rr表示，其中，去除的阵元对应的元素为0，保留的阵元对应的元素为1。

所述应用CLEAN思想清除检测到的目标响应是：将步骤2所计算的所有快拍数据形成波束的幅度平方累加结果P近似为一个快拍数据形成波束的幅度平方累加，首先求得目标T_m对应的幅度

$A\left(T_m\right)=\frac{\mathrm{sqrt}\left(P_{m-1}\left(z_m\right)\right)}{\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{A}_{\mathrm{rr}}\left(n\right)e^{\mathrm{jknd}\sin {\mathrm{\&theta;}}_m}\&CenterDot;e^{-\frac{j2\mathrm{\&pi;}{z}_{m}n}{M}}}$

其中，P_m-1为第m-1次迭代的结果，z_m为P_m-1中最大值对应的波束号，θ_m为步骤(f)

计算的目标T_m的角度，N为阵元数，k＝2π/λ；

在P_m-1清除目标T_m的响应

$P_m\left(q\right)=P_{m-1}\left(q\right)-{|\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}A\left(T_m\right)e^{\mathrm{jknd}\sin {\mathrm{\&theta;}}_m}{A}_{\mathrm{rr}}\left(n\right)\&CenterDot;e^{-\frac{j2\mathrm{\&pi;qn}}{M}}|}^2,q=0\&CenterDot;\&CenterDot;\&CenterDot;M-1.$

有益效果

本发明提出的一种基于CLEAN的稀疏阵列波达方向DOA估计方法，包括以下步骤：(a)利用最小描述长度(MDL)准则来计算目标个数；(b)计算所有快拍数据形成波束的幅度平方累加结果；(c)在计算结果中查找最大值对应的波束号；(d)利用等信号测角法计算“精确波束号”；(e)计算目标角度；(f)判断所有信号是否均已检测，如果是则结束，否则转至(g)；(g)应用CLEAN思想清除检测到的目标响应，转至(c)进行迭代。该方法能够很好的抑制强目标与高旁瓣作用结果对弱目标角度估计的影响，达到很高的测角精度，同时，该方法耗时很短且耗时对阵元数不敏感。

由于本发明的方法利用CLEAN思想能够清除强信号对弱信号检测的影响，并且算法计算比较简单，因此能够在强弱信号共存的情况下，仍然能够快速、精确的估计出稀疏阵列的波达方向。

附图说明

所有附图中FAC表示本发明方法，RM表示Root-MUSIC方法，

图1为示出本发明方法之步骤的流程图；

图2为阵元数都是28的情况下，阵元间距为0.5909λ(λ为波长)的稀疏阵(图中用SA表示)和阵元间距为0.5λ的一维均匀线阵(图中用ULA表示)的方向图对比示意图。

图3为本实例方法与Root-MUSIC方法相对于阵元数的时间对比图。

图4为本实例方法求两个目标(一个SNR为20dB，角度为10°；另一个SNR为40dB，角度为30°)的示例图，其中(a)是最初的响应图，(b)是强目标的响应图，(c)是(a)中去掉(b)得到的响应图。

图5是本实例方法和Root-MUSIC方法相对于SNR的平均均方误差根(RSME)的对比图，目标数为2，目标角度分别为10°和30°，两目标的SNR相等，从10dB变化到70dB，间隔为1dB。

图6是当两目标的SNR相同时(都为20dB)，入射角度比较接近时(其中一个目标角度固定为20°，另一个目标角度从20.15°变化到26°，间隔为0.1°)，本实例方法和Root-MUSIC方法性能对比的示意图，其中，(a)是RMSE对比图，(b)是角度对比图(蓝色实线(angle-1)和虚线(angle-2)表示两个目标的真实角度，红色和绿色分表表示FAC和RM求得的角度值))。

图7是当两目标的SNR不同时时(目标1为20dB，目标2为25dB)，入射角度比较接近时(其中目标1固定为20°，目标2从20.15°变化到26°，间隔为0.1°)，本实例方法和Root-MUSIC方法性能对比的示意图，其中，(a)是RMSE对比图，(b)是角度对比图(蓝色实线(angle-1)和虚线(angle-2)表示两个目标的真实角度，红色和绿色分表表示FAC和RM求得的角度值))

具体实施方式

现结合实施例、附图对本发明作进一步描述：

本发明内容中以一维阵列为例进行描述。

对于一个由N个阵元组成的一维线性阵列天线，在其中随机去除一些阵元可以得到稀疏阵列。为描述方便，用一个包含N个元素的一维数组A_rr表示稀疏阵列，其中，去除的阵元对应的元素为0，保留的阵元对应的元素为1。假定阵元间距为d，有p个入射目标，方向分别为θ_i(i＝1,…p)，各阵元接收到的噪声为零均值高斯白噪声。则接收到得信号可以表示为：

X＝(AS+N)·A_rr (1)式中，X为采样数据，S为目标矢量，N为噪声矢量，A为导向矢量矩阵，具体为：

A＝[a(θ₁)a(θ₂)…a(θ_p)] (2)式中，a(θ_i)为第i个信号的导向矢量，具体为：

$a\left({\mathrm{\&theta;}}_i\right)=[1,e^{\mathrm{jkd}\sin {\mathrm{\&theta;}}_i},\&CenterDot;\&CenterDot;\&CenterDot;,e^{\mathrm{jk}(N-1)d\sin {\mathrm{\&theta;}}_i}{]}^T---\left(3\right)$ 式中k＝2π/λ，d为阵元间距。

现参考图1描述根据本发明的稀疏阵列波达方向估计方法

本发明方法首先需要判断目标的个数。目标个数的准确确定对DOA估计十分重要，如果目标个数选择过大，则会出现虚假目标，如果选择过小，则会漏掉目标。如图1所示，在步骤S101中，本算法中采用最小描述长度(MDL)准则来计算目标个数T。

在形成波束时，对各阵元的采样数据进行加权求和，得到阵列输出(4)

y＝W^HX＝W^H{(AS+N)·A_rr} (4)其中，W＝[W₁,W₂,…W_N]^T为各阵元对应的权值。当W＝a(θ_i)时，波束将指向θ_i方向。

此时，

$y=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\{(\mathrm{AS}+N)\&CenterDot;A_{\mathrm{rr}}\}{e}^{-\mathrm{jknd}\sin {\mathrm{\&theta;}}_i}---\left(5\right)$

若按式(6)进行量化，

$\sin \left({\mathrm{\&theta;}}_i\right)=\frac{\mathrm{\&lambda;q}}{\mathrm{Md}}---\left(6\right)$

则式(5)就成了X(t)的DFT表达式，实际中，经常用FFT来实现。所以对阵列得到的采样值进行M点FFT，就可以得到M个接收波束数据，各波束号与对应的波束指向角度的关系如式(6)所示。

在步骤S102中，判断T是否大于零，如果大于零则进行步骤S103，否则结束计算过程

在步骤S103中，假设第s个快拍数据为x_s，则波束形成结果为：y_s＝FFT(x_s)。对于快拍数为L的数据，得到的所有快拍数据形成波束的幅度平方的累加结果为：

$P=\underset{s=1}{\overset{L}{\mathrm{\&Sigma;}}}{\left|y_s\right|}^2---\left(7\right)$

由FFT形成的两个波束之间的差波束与和波束之比与目标同等信号轴的偏差近似成正比，类似的，先对两波束的幅度求平方，然后再求和波束Σ和差波束Δ，Δ/Σ仍然与目标偏离等信号轴的偏差成正比。因此，利用式(7)求得各波束的累加结果后，P中相邻两波束的差波束与和波束之比与目标偏离等信号轴为线性关系。

在步骤S104中，找出P中最大值对应的波束号p₁，以及该目标对应的次大值的波束号p₂。

$p_2=\left\{\begin{array}{cc}{p}_1-1,& P(p_1-1)>P(p_1+1)\\ p_1+1,& P(p_1-1)\&le;P(p_1+1)\end{array}\right.---\left(8\right)$

在步骤S105中，利用式(9)得到该目标对应的“精确波束号”p_t。

$p_t=\frac{(p_1+p_2)}{2}+k\&CenterDot;\mathrm{sign}\&CenterDot;\min (\frac{|P\left(p_1\right)-P\left(p_2\right)|}{P\left(p_1\right)+P\left(p_2\right)|},0.5)---\left(9\right)$

其中，k为阵列对应的比幅测角斜率，与具体阵列相关，利用式(5)，求得多波束响应y，然后取P＝|y|²，将p₁和p₂选为0…M-1中任意相邻的值，再将p_t取为p₁和p₂之间的值，利用式(9)可求出k。sign为符号位，取值如下：

$\mathrm{sign}=\left\{\begin{array}{cc}-1,& p_1<p_2\\ 1,& p_1\&GreaterEqual;p_2\end{array}\right.---\left(10\right)$

在步骤S106中，利用式(6)可得该目标对应的角度为：

$\mathrm{\&theta;}=\arcsin \left(\frac{p_t\mathrm{\&lambda;}}{\mathrm{Md}}\right)---\left(11\right)$

由于稀疏阵列的旁瓣比较高，产生的假目标在P中对应的值很可能大于真实的弱目标在P中对应的值，所以P中的极大值点有可能对应的是假目标。因此，当MDL求得有T个目标时，在P中取最大的T个极大值点认为是目标是不正确的。

在步骤S107中，为了抑制高旁瓣的影响，应用CLEAN的思想，在每得到一个目标后，在计算结果中清除该目标的响应。精确的做法应该是在每个快拍清除该目标带来的影响，但这样算法复杂度会很大，在实时性要求高的场合达不到***要求。通过仿真实验，我们发现，把所有快拍的计算结果(即P)看作是一个快拍的计算结果，仅作一次清除操作得到的结果与每个快拍做一次清除操作得到的结果十分接近。假设目标T_m在P_m-1(第m-1次迭代的结果)中对应最大波束号为z_m，由式(11)求得的角度为θ_m，则该目标的幅度为：

$A\left(T_m\right)=\frac{\mathrm{sqrt}\left(P_{m-1}\left(z_m\right)\right)}{\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{A}_{\mathrm{rr}}\left(n\right)e^{\mathrm{jknd}\sin {\mathrm{\&theta;}}_m}\&CenterDot;e^{-\frac{j2\mathrm{\&pi;}{z}_{m}n}{M}}}---\left(12\right)$

在P_m-1清除目标T_m的响应后，得到：

$P_m\left(q\right)=P_{m-1}\left(q\right)-{|\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}A\left(T_m\right)e^{\mathrm{jknd}\sin {\mathrm{\&theta;}}_m}{A}_{\mathrm{rr}}\left(n\right)\&CenterDot;e^{-\frac{j2\mathrm{\&pi;qn}}{M}}|}^2,q=0\&CenterDot;\&CenterDot;\&CenterDot;M-1---\left(13\right)$

执行完步骤S107后，转至步骤S102进行迭代。

为说明本发明方法的有效性，下面将本发明方法与Root-MUSIC方法进行仿真实验对比。在本发明后续内容中，FAC表示本发明方法，RM表示Root-MUSIC方法。

不失一般性，我们选用了稀疏阵列A_rr＝[1 1 0 0 0 1 1 0 0 0 0 0 1 0 0 0 1 010 0 1 1 0 1 1 0 0 0 1 1 1 1 0 0 1 0 0 0 1 0 0 0 1 0 1 1 0 0 1 1 1 0 1 1 0 001 1 1]，阵元间距为0.5909λ，所有涉及蒙特卡罗实验的次数(MC)都是200次。实验中所加噪声均为加性高斯白噪声，快拍数为64，所有FFT的点数均为128点。除了实验一外，其他实验都用两个目标进行实验并且加了泰勒窗。对比方法主要是平均均方误差(RMSE)，2个目标对应的计算公式为：

$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{MC}}\left(\underset{i=1}{\overset{\mathrm{MC}}{\mathrm{\&Sigma;}}}\frac{{({\widehat{\mathrm{\&theta;}}}_1-{\mathrm{\&theta;}}_1)}^2+{({\widehat{\mathrm{\&theta;}}}_2-{\mathrm{\&theta;}}_2)}^2}{2}\right)}---\left(14\right)$

实验一：本实验中采用的稀疏阵列有效阵元数(A_rr中为1的个数)为28，将此稀疏阵的方向图与阵元间距为0.5λ的28个均匀线阵的方向图进行对比，结果如图2所示。可以看出，相同(有效)阵元情况下，稀疏阵(对应图中红色)可以获得更窄的主瓣，但均值旁瓣和峰值旁瓣都比较高。

实验二：为对比FAC与RM的时间复杂度，选取阵元数为6至60(间隔为1)的一维均匀线阵，两个目标的入射方向分别为20°和40°，SNR都是20dB。实验结果如图3所示。由图3可以看出在阵元数小于20时，RM耗时较少，在阵元数大于等于20时，FAC耗时较少；同时可以看出，当阵元数增加时，RM的时间快速增加，但FAC的时间增加十分缓慢。可以说明FAC方法的时间对阵元数不敏感，耗时总是很少。

实验三：为验证当强目标旁瓣的响应大于弱目标响应的情况下FAC的性能，本实验中，两目标的SNR分别选为20dB和40dB，入射角度分别为10°和30°。利用式(7)计算的响应图见图4(a)，可以看出最大值点对应的是30°的目标，而10°目标对应的响应并不是一个极值点，因为30°目标带来的旁瓣响应超过了10°目标对应的响应。30°目标对应的响应(也即式(13)中第二项)如图4(b)所示，在图4(a)中去除图4(b)得到图4(c)，可以看出30°目标对应的响应被很好地清除了，10°目标对应的峰值显现了。FAC求得的结果为：10.031°和29.999°，RM计算的结果为：10.044°和29.946°。通过本实验说明，FAC在强目标旁瓣响应大于弱目标响应的情况依然能够精确的估计出所有目标的入射角度。

实验四：两个目标角度选为10°和30°，SNR相等，都是从10dB变化到70dB，间隔为1dB，实验结果如图5所示。由图5可以看出，当两目标入射角度相差较大，且SNR相同时，FAC和RM的精度都很高，FAC优于RM，且FAC与RM都对SNR变化不敏感。

实验五：为验证当两目标SNR相同，角度接近时，算法的性能，两目标的SNR都选为20dB，一个目标的入射角度固定为20°，另一个目标的入射角度从20.15°变化到26°，间隔为0.1°，实验结果如图6所示。由图6(a)可以看出，当SNR相同，两个目标入射角度很接近时，RM和FAC的效果都比较差，在角度十分接近时FAC效果较好，随着角度的增加，逐渐变得RM效果较好，当角度相差较大时，二者效果都很好。

图6(b)中蓝色实线(angle-1)和虚线(angle-2)表示两个目标的真实角度，红色和绿色分别表示FAC和RM求得的角度。可以看出，在两个目标角度相近时，FAC和RM都不是很精确。当两目标角度相差较大时，FAC和RM求得的角度都很精确。

实验六：为验证当两目标SNR不同，角度接近时，算法的性能，第一个目标的SNR取为20dB，第二个目标的SNR取为25dB，第一个目标的入射角度固定为20°，第二个目标的入射角度从20.15°变化到26°，间隔为0.1°，实验结果如图7所示。由图7(a)可以看出，当两个目标的SNR不同时，在大多数情况下FAC的结果优于RM的结果。从图7(b)可以看出，当两目标的SNR不同时，在两个目标角度比较接近时，在大多数情况下，FAC方法几乎可以精确的估计出SNR较高的那个角度，而另一个角度的值和RM估计的基本一致；而RM估计的两个角度都有较大误差。

综合考虑所有仿真实验，可以说明，本发明方法耗时很少且耗时对阵元数不敏感，对目标角度估计的精度很高，能够很好的抑制稀疏阵强目标与高旁瓣作用结果对弱目标角度估计的影响。

Claims (3)

1.一种基于CLEAN的稀疏阵列波达方向DOA估计方法，其特征在于步骤如下：

步骤1：利用最小描述长度MDL准则计算CLEAN的稀疏阵列信号中的目标个数T；

步骤2：T>0时，计算所有快拍数据形成波束的幅度平方累加结果否则结束计算过程；其中：y_s为第s个快拍数据x_s利用FFT形成的波束数据，y_s＝FFT(x_s)；L为快拍数；

步骤3：在幅度平方累加结果P中查找最大值对应的波束号；

步骤4：采用等信号测角法计算“精确波束号”p_t，

$p_t=\frac{(p_1+p_2)}{2}+k^{\&prime;}\&CenterDot;sign\&CenterDot;min(\frac{|P\left(p_1\right)-P\left(p_2\right)|}{|P\left(p_1\right)+P\left(p_2\right)|},0.5)$

其中，p₁为P中最大值对应的波束号，p₂为次大值的波束号，k′为阵列对应的比幅测角斜率，sign为符号位；

$p_2=\left\{\begin{array}{cc}{p}_1-1,& P(p_1-1)>P(p_1+1)\\ p_1+1,& P(p_1-1)\&le;P(p_1+1)\end{array}\right.$

所述

步骤5：计算目标角度，

其中，λ为波长，d为阵元间距，M为FFT点数；

步骤6：应用CLEAN思想清除检测到的目标响应，同时将T值减一，转至步骤2进行迭代。

2.根据权利要求1所述基于CLEAN的稀疏阵列波达方向DOA估计方法，其特征在于：所述稀疏阵列是在均匀间隔阵列中随机去除一些阵元形成的，用数组A_rr表示，其中，去除的阵元对应的元素为0，保留的阵元对应的元素为1。

3.根据权利要求1所述基于CLEAN的稀疏阵列波达方向DOA估计方法，其特征在于：所述应用CLEAN思想清除检测到的目标响应是：将步骤2所计算的所有快拍数据形成波束的幅度平方累加结果P近似为一个快拍数据形成波束的幅度平方累加，首先求得目标T_m对应的幅度

$A\left(T_m\right)=\frac{sqrt\left(P_{m-1}\right(z_m))}{\underset{n=0}{\overset{N-1}{\&Sigma;}}{A}_{rr}\left(n\right)e^{{\mathrm{jkndsin\&theta;}}_m}\&CenterDot;e^{-\frac{j2{\mathrm{\&pi;z}}_{m}n}{M}}}$

其中，P_m-1为第m-1次迭代的结果，z_m为P_m-1中最大值对应的波束号，θ_m为步骤5计算的目标T_m的角度，N为阵元数，k＝2π/λ，A_rr表示稀疏阵列；

在P_m-1清除目标T_m的响应

$P_m\left(q\right)=P_{m-1}\left(q\right)-|\underset{n=0}{\overset{N-1}{\&Sigma;}}A\left(T_m\right)e^{{\mathrm{jkndsin\&theta;}}_m}{A}_{rr}\left(n\right)\&CenterDot;e^{-\frac{j2\mathrm{\&pi;}qn}{M}}{|}^2,q=\mathrm{0...}M-1.$