CN110579737A - 一种杂波环境中基于稀疏阵列的mimo雷达宽带doa计算方法 - Google Patents

Abstract

本发明公开了一种杂波环境中基于稀疏阵列的MIMO雷达宽带DOA计算方法，属于信号处理领域。具体为WCSAB，在该方法中，用Capon波束形成来抑制杂波干扰，然后联合利用不同窄带信号信息，用CS方法估计目标DOA。考虑到DOA估计性能不仅与波束形成权重值有关，还与稀疏阵列结构有关，本发明提出了波束形成权重值和稀疏阵列的联合优化问题，并给出了一种简单的算法求解这一优化问题。本发明提出的方法可以提高杂波环境中目标DOA估计的性能，包括高分辨率和低旁瓣，稀疏阵列降低了***成本和复杂性，以目标DOA估计的贝叶斯均方误差(BMSE)为性能评价指标，通过上述算法设计的稀疏阵列结构与穷举法得到的最优稀疏阵列性能相近，且比嵌套阵和互质阵性能更优。

Description

一种杂波环境中基于稀疏阵列的MIMO雷达宽带DOA计算方法

技术领域

本发明属于信号处理领域，它特别涉及杂波环境中基于稀疏阵列的MIMO雷达宽带DOA 估计问题。

背景技术

MIMO(Multiple Input Multiple Output)雷达是利用多个发射天线同步地发射信号，同时使用多个接收天线接收回波信号，并集中处理的一种新型雷达***。相比于传统的相控阵雷达， MIMO雷达具有明显优势，如更高的分辨率，更好的目标检测、定位和跟踪性能，更好的目标参数估计和识别能力。DOA估计研究是阵列信号处理中的一项重要内容，其应用涉及雷达、通信、声呐、射电天文、勘测、地震以及生物医学等领域。目前已有多种经典的DOA估计方法，如多信号分类(Multiple Signal Classification,MUSIC)，基于旋转不变技术的信号参数估计(Estimation of Signal Parameters via Rotational invariancetechniques,ESPRIT)等。近年来，压缩感知(Compressive Sensing,CS)理论得到了国内外学者的广泛关注，相对于传统的方法，基于CS的MIMO雷达DOA估计在较少采样数据和低信噪比条件下拥有更好的估计性能。

根据传统的阵列信号处理理论，为了保证DOA估计的唯一性，阵列中相邻阵元间距要小于或等于入射信号半波长，满足这一条件的阵列称为满阵。阵列的空间分辨率与阵列孔径有关，提高分辨率需要增大阵列孔径，在满阵中，这意味着需要更多的天线。然而，由于实际中软硬件资源的约束，天线的数目通常是受限的。为了增大阵列孔径，同时不增加天线个数，稀疏阵列引起了广泛的关注。当目标在观测空间是稀疏时，稀疏阵列能够准确的估计出目标DOA。然而，在杂波环境中，目标在观测空间的稀疏性会遭到破坏，从而导致DOA估计性能下降。

考虑到宽带信号具有信息量大、抗干扰能力强、分辨率高等优点，比较有代表性的宽带 DOA估计方法有非相干信号子空间(ISSM)和相干信号子空间(CSSM)等。ISSM将宽带信号分成频带上的多个窄带信号，然后分别对每个窄带信号进行处理，最后对所有窄带信号处理的结果求平均得到最终的估计结果，这种方法计算量大，且在低信噪比时估计性能差。CSSM 是通过聚焦将不同频率窄带信号的协方差矩阵变换到一个参考频率，然后再用窄带估计方法得到最终结果。但是CSSM需要对目标DOA预估，且性能受预估精度的影响很大。

发明内容

本发明提供了一种杂波环境中基于稀疏阵列的MIMO雷达宽带DOA估计方法，具体为 WCSAB(wideband compressive sensing after beamforming)，在该方法中，用Capon波束形成来抑制杂波干扰，然后联合利用不同窄带信号信息，用CS方法估计目标DOA。考虑到DOA 估计性能不仅与波束形成权重值有关，还与稀疏阵列结构有关，本发明提出了波束形成权重值和稀疏阵列的联合优化问题，并给出了一种简单的算法求解这一优化问题。

本发明技术方案为一种杂波环境中基于稀疏阵列的MIMO雷达宽带DOA计算方法，该方法包括：

步骤1：设发射天线位置确定，放置接收天线的可行域为[0,D_r]，为了简化分析，将可行域以间隔Δ_r离散化为N_r个格点，且有N个接收天线放置在其中一些格点上，N＜＜N_r；

步骤2：建立MIMO雷达回波信号模型，得到回波信号时域采样数据 n＝1,...,N_r和p＝1,...,L，其中p表示时域快拍，L为快拍数；

步骤3：对接收信号进行L点离散傅里叶变换得到频域数据，即

并将N_r个格点的数据表示成矢量形式，即y[l]＝[y₁[l],...,y_Nr[l]]^T，其中p＝1,...,L，l＝1,...,L；

步骤4：将目标角度观测区域离散化为G个格点θ₁,...,θ_G，K＜＜G，其中K表示目标个数，将信号模型表示成稀疏形式：

y[l]＝Φ[l]x+c[l]+u[l]

其中这里a_r(θ,f_l)表示接收导向矢量， a_t(θ,f_l)表示发射导向矢量，s[l]表示频域发射信号，x＝[x₁,...,x_G]^T是K稀疏的，也就是x只有 K个非零元素，且非零元素的值和位置为目标反射系数和DOA，c[l]表示杂波，u[l]表示噪声；

步骤5：将波束形成权重矢量w_g,l作用到y[l]上得到波束形成输出结果：

将r_g,l(g＝1,...,G，l＝1,...,L)表示成矢量：

r＝[r_1,1,...,r_G,1,...,r_1,L,...,r_G,L]^T

＝W_rΦx+W_rc+W_ru

其中权值矩阵W_r＝Diag{W₁,...,W_l,...,W_L}是一个块对角矩阵；

且有Φ＝[Φ^T[1],...,Φ^T[L]]^T，c＝[c^T[1],...,c^T[L]]^T表示杂波， u＝[u^T[1],...,u^T[L]]^T表示噪声；

步骤6：基于CS理论，通过基寻踪去燥来重构稀疏向量x；

其中η≥0是正则化参数；

步骤7：对步骤6得到的解中的元素值按照从大到小排序，排序后各个元素相应的格点表示为{θ₍₁₎,...,θ_(G)}，那么DOA估计结果可以表示为；

步骤8：基于最小化贝叶斯均方误差求解最优W_r，建立以下优化问题

s.t.W_r＝Diag{W₁,...,W_L}

||W||₀＝N

W＝[w_1,1,...,w_1,L,...,w_G,1,...,w_G,L]

其中，真实目标的DOA矢量θ_T是随机的，表示对θ_T求期望，表示θ_T确定时，DOA估计的均方误差，w_g,l表示权重矢量，其中g＝1,...,G，l＝1,...,L；

步骤9：优化求解步骤8提出的问题，得到最优的W_r。

进一步的，所述步骤9的具体方法为：

步骤1：初始化：迭代次数j＝1，根据公式计算波束形成权重值其中R_c(f_l)为杂波c[l]的协方差矩阵；在每次迭代过程中，随机产生一组格点选择矢量{z₁,...,z_α}，对于给定的z；

步骤2：重复步骤3到步骤6的迭代过程：

步骤3：随机产生一组格点选择矢量{z₁,...,z_α}；

步骤4：根据公式w_g,l＝z⊙ξ_g,l计算并构成

根据公式计算r_g,l，将r_g,l表示成矢量r；

将和r代入公式得到x重构结果和目标DOA估计结果

根据公式得到BMSE

步骤5：基于最小BMSE得到

步骤6：根据基于更新中相应的权重值得到并令j＝j+1；其中为Dc[l]的协方差矩阵，

步骤7：当时，迭代停止，输出最优天线选择e₀为事先设定的阈值。

本发明提出的方法可以提高杂波环境中目标DOA估计的性能，包括高分辨率和低旁瓣，稀疏阵列降低了***成本和复杂性，以目标DOA估计的贝叶斯均方误差(BMSE)为性能评价指标，通过上述算法设计的稀疏阵列结构与穷举法得到的最优稀疏阵列性能相近，且比嵌套阵和互质阵性能更优。

附图说明

图1给出了所有可能稀疏阵列结构下的BMSE以升序排列的结果，为了对比，图1也给出了嵌套阵(nested array)和互质阵(co-prime array)的结果。

图2(a)给出了基于最小BMSE条件下最优的稀疏阵列结构，图2(b)给出了根据本发明所提出的算法得到的稀疏阵列结构。

图3给出了使用WCSAB方法时，不同稀疏阵列结构的DOA估计结果。

图4考虑单目标情况，分别用WCSAB和WCT(wideband Capon technique)方法时，不同阵列结构的DOA估计结果。

图5为考虑双目标情况，分别用WCSAB和WCT(wideband Capon technique)方法时，不同阵列结构的DOA估计结果。

具体实施方式

为了方便描述，首先进行如下定义：

黑体大写字母表示矩阵，黑体小写字母表示矢量，(·)^*为共轭，(·)^T为转置，(·)^H为共轭转置，||x||₀和||x||₁分别表示向量x的l₀范数和l₁范数，||W||₀表示矩阵W非零行的个数，Diag{·}表示块对角矩阵，diag_r{·}表示去掉零行之后的对角矩阵，表示相对于θ的期望，I_N为N阶的单位阵，1为全1矢量，符号⊙表示哈达玛积。

考虑一个共置MIMO雷达***，发射天线和接收天线都放置在二维笛卡尔坐标系的横轴上。假设有M个发射天线，且在横轴上的位置已知，为d_t,m(m＝1,...,M)。假设放置接收天线的可行域为[0,D_r]，为了简化分析，将可行域以间隔Δ_r离散化为N_r个格点，接收天线放置在这些格点上。由于天线个数的约束，假设雷达***只有N(N＜＜N_r)个可用的接收天线。令表示第m个发射天线发射的宽带信号，频率范围为[-B_m/2,B_m/2]，其中p 表示时域快拍，T_s表示采样周期，L表示快拍数。假设K个远场点目标的DOA为θ_T,k(k＝1,...,K)，那么在第n个格点接收到的信号为

其中f_c表示载频，β_k表示第k个目标的反射系数，且假设是确定未知的。令第一个发射天线和第一个格点作为参考，那么τ_Tt,m,k＝(d_t,m-d_t,1)sinθ_T,k/c表示信号从第m个发射天线到第k个目标时，相对于参考阵元的时延，τ_Tr,n,k＝(n-1)Δ_rsinθ_T,k表示信号从第k个目标到第n个格点时，相对于第一个格点的时延。Q表示杂波散射体的个数，γ_q(q＝1,...,Q)表示杂波散射体的反射系数，并且假设它们之间是独立同分布(iid)的高斯随机变量。τ_Ct,m,q＝(d_t,m-d_t,1)sinθ_C,q/c 表示信号从第m个发射天线到第q个杂波散射体时，相对于参考阵元的时延，τ_Cr,n,q＝(n-1)Δ_rsinθ_C,q表示信号从第q个杂波散射体到第n个格点时，相对于第一个格点的时延，θ_C,q表示第q个杂波散射体相对于阵列的方向。是方差为σ²的高斯白噪声。

通过对时域离散信号进行L点离散傅里叶变换(DFT)，可得到在频率点f_l＝lf_s(l＝1,...,L) 的频域数据，其中f_s为频率采样间隔，f_l∈[-B/2,B/2]且信号在频率f_l处的 DFT结果为

其中s_m[l]和u_n[l]分别表示发射信号和噪声的DFT。令和分别表示在角度θ、频率f_l处的接收导向矢量和发射导向矢量。将N_r个格点接收到的信号表示成矢量

其中

在CS框架下，为了估计K个目标的DOA θ_T,k(k＝1,...,K)，将目标角度观测区域离散化为G(K＜＜G)个格点θ₁,...,θ_G，假设离散误差可以忽略，即目标正好落在格点上。那么(3)式可以表示为

y[l]＝Φ[l]x+c[l]+u[l] (4)

其中矢量x＝[x₁,...,x_G]^T是K稀疏的，也就是x只有K个非零元素，且非零元素的值和位置为目标反射系数和DOA，可以表示为

CS理论利用x的稀疏性来估计目标DOA，然而，这种稀疏性在杂波环境中会遭到破坏，从而降低DOA估计的性能。为了抑制杂波的干扰，在接收端采用波束形成的方法。令表示在方向θ_g、频率f_l处的波束形成权重矢量，且非零元素的位置表示选择放置天线的格点。由于只有N个可用的接收天线，因此要求权重矢量满足||w_g,l||₀＝N。波束形成的输出由下式给出

将r_g,l(g＝1,...,G和l＝1,...,L)表示成一个GL×1的矢量

其中Φ＝[Φ^T[1],...,Φ^T[L]]^T，c＝[c^T[1],...,c^T[L]]^T，u＝[u^T[1],...,u^T[L]]^T，W_r＝Diag{W₁,...,W_L}，根据(7)式，DOA估计问题可以转化为稀疏信号重构问题，基于CS理论，K稀疏的矢量x可以通过基寻踪去燥(BPDN)来重构

其中η≥0是正则化参数，对于(8)式这个优化问题，可以使用CVX工具包求解。令表示上式的解，为目标DOA的估计结果。考虑目标DOA矢量θ_T＝[θ_T,1,...,θ_T,K]^T是随机的情况，那么平均估计性能可以由贝叶斯均方误差(BMSE)给出

由(9)式可知，DOA估计的性能与矩阵W_r有关，为了使性能最优，给出以下优化问题

(10)式中最后两个约束是为了保证对于不同的g和l，w_g,l中非零元素的位置是相同的。由于 w_g,l中非零元素的位置表示选中相应的格点放置天线，因此(10)式是一个权重值和稀疏阵列结构的联合优化问题。

考虑到(10)式是一个NP-hard问题，提出一种简单的算法来求解该优化问题。该算法首先在权值给定时优化稀疏阵列结构，然后再去更新权值用于下一次迭代。首先解释如何在权值给定时优化稀疏阵列结构。首先，定义一个格点选择矢量其中z_n∈{0,1}，只有元素为1时表示选择相应的格点放置天线，由于只有N个可用接收天线，要求||z||₀＝N。第一次迭代的权值由满阵情况Capon波束形成给出

其中R_c(f_l)为杂波c[l]的协方差矩阵。在每次迭代过程中，随机产生一组格点选择矢量{z₁,...,z_α}，对于给定的z，有

w_g,l＝z⊙ξ_g,l (12)

对于不同的z，可以得到不同的w_g,l和W_r(W_r由w_g,l构成)，由(9)式可知，BMSE与W_r有关，因此可以发现，BMSE与z也相关，表示为e(z)。基于最小BMSE，可以得到最优的格点选择矢量z_op。

下面基于z_op更新ξ_g,l中相应的权重值，令则ξ_g,l中由z_op选中的元素值通过下式更新

其中为Dc[l]的协方差矩阵，当BMSE e(z_op)小于某个阈值e₀时，迭代停止。详细的算法由表1给出。

表1求解优化问题的迭代算法

为了抑制杂波的干扰，在接收端采用波束形成的方法。表示在方向θ_g、频率f_l处的波束形成权重矢量，且非零元素的位置表示选择放置天线的格点，由于只有N个可用的接收天线，因此要求权重矢量满足||w_g,l||₀＝N。对于不同的g和l，天线的阵列结构应是相同的，也就是w_g,l中非零元素的位置是相同的。为了表示这一约束，构造矩阵

W＝[w_1,1,...,w_1,L,...,w_G,1,...,w_G,L] (14)

且满足||W||₀＝N，即矩阵非零行的个数为N，通过这个约束，即可满足对于不同的g和l，w_g,l中非零元素的位置相同。将权重矢量w_g,l作用于接收信号y[l]，根据(6)式得到波束形成输出 r_g,l，将其表示成G×1的矢量

由上式可以发现，对于不同的频率f_l，向量x稀疏性一样。为了联合利用不同频率的信号信息，将r_l(l＝1,...,L)表示成GL×1的矢量，即(7)式。通过将DOA估计问题可以转化为稀疏信号重构问题，根据(8)式，可以得到稀疏向量x的重构结果则中最大的K个元素的位置为目标DOA的估计结果，表示为对中的元素值按照从大到小排序，排序后各个元素相应的格点表示为{θ₍₁₎,...,θ_(G)}，那么DOA估计结果可以表示为

满阵波束形成权重矢量ξ_g,l可以使方向θ_g、频率f_l处的信号无失真通过，同时抑制其他方向的干扰和噪声，表示为

上式的最优解即为(11)式。

关于杂波环境中基于稀疏阵列的MIMO雷达宽带DOA估计，给出了两个仿真实例，参数设置如下：假设D_r＝11λ/2，其中λ表示信号最高频率对应的波长。将可放置接收天线的可行域以Δ_r＝λ/2为间隔离散为12个格点。假设MIMO雷达***可用的发射和接收天线个数为M＝N＝6，且发射端阵列结构确定已知。

为了简化分析，假设发射信号带宽相同，即B_m＝200MHz(m＝1,...,M)，载频为1GHz。

将目标角度观测域离散为41个格点-20°,-19°,...,20°。

杂波由250个散射体构成，分布的角度为-90°,-90°+180°/250,...,90°。

定义信噪比和信杂比不失一般性，假设目标反射系数为1，SNR和SCR分别设置为-5dB和-30dB。

在仿真1中，假设目标在角度观测域离散化后的格点上是均匀随机分布的。为了保证阵列孔径不改变，令那么共有种不同的稀疏阵列结构。图1给出了所有可能稀疏阵列结构下的BMSE以升序排列的结果，菱形表示最小BMSE条件下最优的稀疏阵列结构，正方形表示根据表1所给出的算法得到的稀疏阵列结构，可以看到这两种结构的性能相近，这两种稀疏阵列的具体结构由图2给出。为了对比，图1也给出了嵌套阵(nested array) 和互质阵(co-prime array)的结果，分别由星号和圆形标出，可以看到有多种稀疏阵列结构的性能都比嵌套阵和互质阵更好。假设只有一个目标，DOA为-14°，图3给出了使用WCSAB 方法时，上述四种稀疏阵列结构的DOA估计结果。从图中可以看到，最优稀疏阵列和由算法得到的稀疏阵列可以准确估计目标DOA，而嵌套阵和互质阵估计有误差。

在仿真2中，比较了WCSAB与WCT(wideband Capon technique)这两种宽带DOA估计方法的性能。WCT方法属于ISSM的一种，它对每一个窄带信号用Capon方法得到相应的DOA估计结果，然后对所有的结果求平均得到最终的估计结果。图4考虑单目标情况，目标DOA为10°，图5考虑双目标情况，目标DOA为6°和10°。图中实线表示以λ/2为间隔的满阵结构，虚线表示最优稀疏阵列结构，点划线表示根据算法得到的稀疏阵列结构。从图4和图5可以看出，相比于满阵结构，两种稀疏阵列所对应的主瓣宽度更窄，也就是说稀疏阵列的分辨率更高。但是稀疏阵列会导致较高的旁瓣，并且这一问题在双目标情况下更严重。通过比较，可以发现， WCSAB方法的旁瓣比WCT低。由图5也可以看到，WCSAB方法在两种稀疏阵列中都可以准确的估计出目标DOA，而WCT估计有误差。通过对比，WCSAB方法的性能更好。

Claims (2)

1.一种杂波环境中基于稀疏阵列的MIMO雷达宽带DOA计算方法，该方法包括：

步骤1：设发射天线位置确定，放置接收天线的可行域为[0,D_r]，为了简化分析，将可行域以间隔Δ_r离散化为N_r个格点，且有N个接收天线放置在其中一些格点上，N＜＜N_r；

步骤2：建立MIMO雷达回波信号模型，得到回波信号时域采样数据n＝1,...,N_r和p＝1,...,L，其中p表示时域快拍，L为快拍数；

步骤3：对接收信号进行L点离散傅里叶变换得到频域数据，即

并将N_r个格点的数据表示成矢量形式，即其中p＝1,...,L，l＝1,...,L；

步骤4：将目标角度观测区域离散化为G个格点θ₁,...,θ_G，K＜＜G，其中K表示目标个数，将信号模型表示成稀疏形式：

y[l]＝Φ[l]x+c[l]+u[l]

其中这里a_r(θ,f_l)表示接收导向矢量，a_t(θ,f_l)表示发射导向矢量，s[l]表示频域发射信号，x＝[x₁,...,x_G]^T是K稀疏的，也就是x只有K个非零元素，且非零元素的值和位置为目标反射系数和DOA，c[l]表示杂波，u[l]表示噪声；

步骤5：将波束形成权重矢量w_g,l作用到y[l]上得到波束形成输出结果：

将r_g,l(g＝1,...,G，l＝1,...,L)表示成矢量：

r＝[r_1,1,...,r_G,1,...,r_1,L,...,r_G,L]^T

＝W_rΦx+W_rc+W_ru

其中权值矩阵W_r＝Diag{W₁,...,W_l,...,W_L}是一个块对角矩阵；

且有Φ＝[Φ^T[1],...,Φ^T[L]]^T，c＝[c^T[1],...,c^T[L]]^T表示杂波，u＝[u^T[1],...,u^T[L]]^T表示噪声；

步骤6：基于CS理论，通过基寻踪去燥来重构稀疏向量x；

其中η≥0是正则化参数；

步骤7：对步骤6得到的解中的元素值按照从大到小排序，排序后各个元素相应的格点表示为{θ₍₁₎,...,θ_(G)}，那么DOA估计结果可以表示为；

步骤8：基于最小化贝叶斯均方误差求解最优W_r，建立以下优化问题

s.t.W_r＝Diag{W₁,...,W_L}

||W||₀＝N

W＝[w_1,1,...,w_1,L,...,w_G,1,...,w_G,L]

其中，真实目标的DOA矢量θ_T是随机的，表示对θ_T求期望，表示θ_T确定时，DOA估计的均方误差，w_g,l表示权重矢量，其中g＝1,...,G，l＝1,...,L；

步骤9：优化求解步骤8提出的问题，得到最优的W_r。

2.如权利要求1所述的一种杂波环境中基于稀疏阵列的MIMO雷达宽带DOA计算方法，其特征在于所述步骤9的具体方法为：

步骤9.1：初始化：迭代次数j＝1，根据公式计算波束形成权重值其中R_c(f_l)为杂波c[l]的协方差矩阵；在每次迭代过程中，随机产生一组格点选择矢量{z₁,...,z_α}，对于给定的z；

步骤9.2：重复步骤9.3到步骤9.6的迭代过程：

步骤9.3：随机产生一组格点选择矢量{z₁,...,z_α}；

步骤9.4：根据公式w_g,l＝z⊙ξ_g,l计算并构成

根据公式计算r_g,l，将r_g,l表示成矢量r；

将和r代入公式得到x重构结果和目标DOA估计结果

根据公式得到

步骤9.5：基于最小BMSE得到

步骤9.6：根据基于更新中相应的权重值得到并令j＝j+1；其中 为Dc[l]的协方差矩阵，

步骤9.7：当时，迭代停止，输出最优天线选择e₀为事先设定的阈值。