CN103760546B - 一种雷达用低空目标波达方向估计方法 - Google Patents

Abstract

本发明属于雷达波达方向估计技术领域，公开了一种雷达用低空目标波达方向估计方法。该雷达用目标波达方向估计方法，包括以下步骤：雷达的天线阵列接收回波信号；对所述回波信号进行采样，得到数据矢量矩阵；对数据矢量矩阵进行奇异值分解，生成对应的降维数据矢量、以及噪声子空间矩阵；得出导向矢量阵；根据噪声子空间矩阵、以及导向矢量阵，得出对应的多径衰减矢量；构造复合导向矢量；构造感知矩阵；得出正则化参数；通过求解对应的数学模型得出波达方向。

Description

一种雷达用低空目标波达方向估计方法

技术领域

本发明属于雷达波达方向估计技术领域，特别涉及一种雷达用低空目标波达方向估计方法。

背景技术

雷达用于远程预警时，作用距离较远，由于远距低空目标的飞行高度多在几千米以下的低空领域，强烈的地海杂波将导致回波信号的信杂比降低，从而难以准确得到波达方向（DOA）。同时，雷达在探测和跟踪低空目标时，不仅接收到直达波信号，还会接收到多径反射回波信号，导致多径效应。由于直达波信号和反射信号高度相关，且两回波入射角度大多相差小于一个半功率的波束宽度，因此，二者为相干信号源。传统的单脉冲测角方法和波束扫描法难以区分出目标真实回波信号，无法得到有效的波达方向。因此，能够分辨相干信号的高分辨算法成为该领域研究的热点。

目前，许多学者提出了一些方法来解决低空目标的波达方向检测问题。最大似然估计算法和广义MUSIC算法虽然对相干信号具有高分辨力，但都需要二维空间谱搜索，计算量大；基于多快拍下稀疏矩阵重构的L1-SVD角度估计算法和矩阵束解相干算法虽然降低了运算量，但这两种算法只适用在高信噪比的条件下，在低信杂噪比下其分辨能力仍然较差。

发明内容

本发明的目的在于提出一种雷达用低空目标波达方向估计方法。该雷达用低空目标波达方向估计方法可应用于目标定向和雷达追踪，运算量小，并具有易于工程实现的特点。

为实现上述技术目的，本发明采用如下技术方案予以实现。

一种雷达用低空目标波达方向估计方法，包括以下步骤：

S1：雷达的天线阵列接收回波信号；对所述回波信号进行采样，得到数据矢量矩阵X；对数据矢量矩阵X进行奇异值分解，根据分解结果，生成对应的降维数据矢量X_c、以及噪声子空间矩阵E_n；

S2：根据雷达搜索角度范围的下限、雷达搜索角度范围的上限、以及雷达测角精度，得出角度采样点个数M、以及对应的第1个采样角度至第M个采样角度，计算与第i目标采样角度θ_1i对应的第i镜像采样角度θ_2i，i取1至M；根据第i目标采样角度θ_1i和第i镜像采样角度θ_2i，得出对应的导向矢量阵A(θ_1i,θ_2i)；

S3：根据噪声子空间矩阵E_n、以及导向矢量阵A(θ_1i,θ_2i)，得出对应的多径衰减矢量构造复合导向矢量b(θ_1i)，构造感知矩阵B，B=[b(θ₁₁),...,b(θ_1M)]；

S4：得出满足以下关系式的正则化参数μ：μ²≥Mσ²，σ²为所述回波信号中的噪声信号的平均功率；建立如下范数约束式：

$\underset{{\mathrm{\ψ}}_c}{\min }{\left|\right|{\mathrm{\ψ}}_c\left|\right|}_1$

$s.t.{\left|\right|X_c-B{\mathrm{\ψ}}_c\left|\right|}_2^2\≤{\mathrm{\μ}}^2$

其中，Ψ_c为稀疏矢量，||·||₁表示L1范数，||·||₂表示L2范数，为：||Ψ_c||₁最小时的Ψ_c；然后，针对所述范数约束式，求解出Ψ_c；根据以下公式计算波达方向其中，Ψ_c(i)为求解出的Ψ_c中的第i个元素，abs(Ψ_c(i))表示Ψ_c(i)的模值，表示min(abs(Ψ_c(i)))对应的θ_1i。

本发明的特点和进一步改进在于：

在步骤S1中，所述雷达的天线阵列包括N个阵元，每个阵元在接收到回波信号之后，进行L次采样；第n个阵元在进行第l次采样后得到的数据为x_n(l)，其中，n取1至N，l取1至L；然后构造数据矢量矩阵X，数据矢量矩阵X的行数为N，列数为L，数据矢量矩阵X的第n行第l列的元素为x_n(l)；

在步骤S1中，按照以下公式对数据矢量矩阵X进行奇异值分解：X=UΛV^H，其中，U为N×N阶酉矩阵，Λ为N×L阶矩阵，V为L×L阶酉矩阵，V^H为V的共轭转置；然后根据以下公式生成对应的降维数据矢量X_c：

X_c=UΛY=XVY

其中，Y为N×1维向量，Y中第一行第一列元素为1，其余元素均为0；在N阶酉矩阵U中，取出后N-1列元素构成噪声子空间矩阵E_n。

在步骤S2中，根据以下公式计算得出角度采样点个数M：其中，α为雷达搜索角度范围的下限，β为雷达搜索角度范围的上限，Δθ为雷达测角精度；然后根据以下公式得出第i目标采样角度θ_1i：

θ_1i=α+(i-1)·Δθ；

然后，根据以下公式计算与第i镜像采样角度θ_2i：

${\mathrm{\θ}}_{2i}=-\arctan \frac{R\sin {\mathrm{\θ}}_{1i}+2(h_a-h_s)}{R\cos {\mathrm{\θ}}_{1i}},$

其中，R为目标与雷达之间的距离，h_a为雷达高度，h_s为多径效应的反射面高度，本发明实施例中，h_s等于0。然后，构造目标导向列向量a(θ_1i)和镜像导向列向量a(θ_2i)，

$a\left({\mathrm{\θ}}_{1i}\right)=\left[\begin{array}{c}1\\ e^{-j\frac{2\mathrm{\πf}}{c}d\sin {\mathrm{\θ}}_{1i}}\\ .\\ .\\ .\\ e^{-j\frac{2\mathrm{\πf}}{c}(N-1)d\sin {\mathrm{\θ}}_{1i}}\end{array}\right]a\left({\mathrm{\θ}}_{2i}\right)=\left[\begin{array}{c}1\\ e^{-j\frac{2\mathrm{\πf}}{c}d\sin {\mathrm{\θ}}_{2i}}\\ .\\ .\\ .\\ e^{-j\frac{2\mathrm{\πf}}{c}(N-1)d\sin {\mathrm{\θ}}_{2i}}\end{array}\right]$

其中，f为雷达频率，d为雷达阵元间距，c为光速，则导向矢量阵A(θ_1i,θ_2i)为：A(θ_1i,θ_2i)=[a(θ_1i)a(θ_2i)]。

在步骤S3中，按照以下公式计算噪声投影矩阵P：表示E_n的共轭转置；按照以下公式得出多径衰减矢量

$\hat{p}({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})=\frac{{\left[A^H({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})\mathrm{PA}({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})\right]}^{-1}w}{w^H{\left[A^H({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})\mathrm{PA}({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})\right]}^{-1}w}$

其中，w=[1,0]^T。

在步骤S4中，针对所述范数约束式，利用凸优化设计工具箱CVX，求解出Ψ_c。

本发明的有益效果为：本发明充分利用了多径衰减系数和角度关系，在有效提高角度估计精度的同时，减小角度搜索范围地。通过将复合导向矢量利用到稀疏矩阵重构中，提高了信息矢量的稀疏性，在低信噪比下对多径信号具有较高的角度分辨能力。

附图说明

图1为本发明的一种雷达用低空目标波达方向估计方法的流程示意图；

图2为本发明和L1-SVD算法在波达方向均方根误差上的对比示意图；

图3为利用本发明检测波达方向时的多径衰减系数估计误差随信噪比变化示意图。

具体实施方式

下面结合附图对本发明作进一步说明：

参照图1，为本发明的一种雷达用低空目标波达方向估计方法的流程示意图。该雷达用低空目标波达方向估计方法包括以下步骤：

S1：雷达的天线阵列为N个等距离的阵元（即N元线阵），N为大于1的自然数，N个等距离的阵元包括第1阵元至第N阵元；每个阵元接收目标的回波信号，然后对回波信号进行第1次采样至第L次采样（数字采样），L为大于1的自然数。第n个阵元在进行第l次采样后得到的数据为x_n(l)，其中，n取1至N，l取1至L。

根据x_n(l)，获得数据矢量矩阵X（N×L阶矩阵），数据矢量矩阵X的第n行第l列的元素为x_n(l)，也就是说数据矢量矩阵X为：

然后，按照以下公式对数据矢量矩阵X进行奇异值分解：X=UΛV^H，其中，U为N×N阶酉矩阵，Λ为N×L阶矩阵（X的奇异值从大到小依次排列在Λ的对角线上），V为L×L阶酉矩阵，V^H为V的共轭转置。

根据N阶酉矩阵U和N×L阶矩阵Λ，计算出N×1维的降维数据矢量X_c：

X_c=UΛY=XVY

其中，Y为N×1维向量，Y中第一行第一列元素为1，其余元素均为0。

在N×N阶酉矩阵U中，取出后N-1列元素构成噪声子空间矩阵E_n。

S2：根据以下公式计算得出角度采样点个数M：其中，α为雷达搜索角度范围的下限，β为雷达搜索角度范围的上限，Δθ为雷达测角精度。

根据以下公式得出第i目标采样角度θ_1i：

θ_1i=α+(i-1)·Δθ；

根据以下公式计算与第i镜像采样角度θ_2i：

${\mathrm{\θ}}_{2i}=-\arctan \frac{R\sin {\mathrm{\θ}}_{1i}+2(h_a-h_s)}{R\cos {\mathrm{\θ}}_{1i}},$

其中，R为目标与雷达之间的距离，h_a为雷达高度，h_s为回波反射面高度（h_s等于0）；目标与雷达之间的距离R的获取过程如下，雷达利用脉冲压缩技术发射信号，并利用回波信号计算出目标与雷达之间的距离R。

构造目标导向列向量a(θ_1i)和镜像导向列向量a(θ_2i)：

$a\left({\mathrm{\θ}}_{1i}\right)=\left[\begin{array}{c}1\\ e^{-j\frac{2\mathrm{\πf}}{c}d\sin {\mathrm{\θ}}_{1i}}\\ .\\ .\\ .\\ e^{-j\frac{2\mathrm{\πf}}{c}(N-1)d\sin {\mathrm{\θ}}_{1i}}\end{array}\right]a\left({\mathrm{\θ}}_{2i}\right)=\left[\begin{array}{c}1\\ e^{-j\frac{2\mathrm{\πf}}{c}d\sin {\mathrm{\θ}}_{2i}}\\ .\\ .\\ .\\ e^{-j\frac{2\mathrm{\πf}}{c}(N-1)d\sin {\mathrm{\θ}}_{2i}}\end{array}\right]$

其中，j表示虚部，f为雷达频率，d为雷达阵元间距，c为光速。

构造导向矢量阵A(θ_1i,θ_2i)为：A(θ_1i,θ_2i)=[a(θ_1i)a(θ_2i)]，由此看出，A(θ_1i,θ_2i)为N×2阶矩阵。

S3：按照以下公式计算噪声投影矩阵P：表示E_n的共轭转置。根据噪声子空间矩阵E_n、以及导向矢量阵A(θ_1i,θ_2i)，得出对应的多径衰减矢量 $\hat{p}({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i}):$

$\hat{p}({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})=\frac{{\left[A^H({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})\mathrm{PA}({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})\right]}^{-1}w}{w^H{\left[A^H({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})\mathrm{PA}({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})\right]}^{-1}w}$

其中，w=[1,0]^T。

构造复合导向矢量b(θ_1i)， $b\left({\mathrm{\θ}}_{1i}\right)=A({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})\hat{p}({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i}).$

构造感知矩阵B，B=[b(θ₁₁),...,b(θ_1M)]。

S4：得出满足以下关系式的正则化参数μ：μ²≥Mσ²，σ²为所述回波信号中的噪声信号在时域中的平均功率，以下举例说明σ²的获得方式：将目标的回波信号进行离散傅里叶变换，使空中目标的回波信号转变为频域信号，根据频域信号，估计噪声信号在频域中的平均功率；然后将所述噪声信号在频域中的平均功率除以离散傅里叶变换点数，得出噪声信号在时域中的平均功率σ²。

然后建立如下范数约束式：

$\underset{{\mathrm{\ψ}}_c}{\min }{\left|\right|{\mathrm{\ψ}}_c\left|\right|}_1$

$s.t.{\left|\right|X_c-B{\mathrm{\ψ}}_c\left|\right|}_2^2\≤{\mathrm{\μ}}^2$

其中，Ψ_c为稀疏矢量，||·||₁表示L1范数，||·||₂表示L2范数，为：||Ψ_c||₁最小时的Ψ_c；然后，针对所述范数约束式，利用凸优化设计工具箱CVX，求解出Ψ_c；根据以下公式计算波达方向其中，Ψ_c(i)为求解出的Ψ_c中的第i个元素（从上至下），abs(Ψ_c(i))表示Ψ_c(i)的模值，表示min(abs(Ψ_c(i)))对应的θ_1i。至此，得到了雷达的波达方向（即目标角度）

下面通过一个仿真实施例对本发明作进一步说明：

1)仿真条件说明：

为不失一般性，假设雷达的天线阵列为一个俯仰方向上的均匀线阵，阵元数为20，阵元间距为半波长，雷达载频为10GHz，雷达高度为10m，接收到的信号均为窄带信号，噪声为与信号不相关的高斯白噪声，快拍数（采样数）为50，多径效应的反射面高度为0，目标高度为100m，目标与雷达距离为6km，进行100次蒙特卡洛实验。根据多径模型通过计算可得此时目标直射回波入射角约为15mrad，镜像反射回波入射角约为-18.3mrad。波达方向的均方根误差采用以下公式计算：

$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{Monte}}\underset{t=1}{\overset{\mathrm{Monte}}{\mathrm{\Σ}}}{({\widehat{\mathrm{\θ}}}_t-{\mathrm{\θ}}_{\mathrm{real}})}^2},$

其中，Monte为蒙特卡洛实验次数，因此Monte=100，θ_real为目标真实角度，表示第t次蒙特卡洛实验角度估计值；多径衰减系数估计误差采用如下公式计算：

${\mathrm{\β}}_{\mathrm{err}}=\frac{1}{\mathrm{Monte}}\left(\underset{t=1}{\overset{\mathrm{Monte}}{\mathrm{\Σ}}}\mathrm{abs}({\widehat{\mathrm{\β}}}_t-\mathrm{\β})\right)$

其中，β为多径衰减系数的真实值，表示第i次蒙特卡洛实验系数估计值。

2)仿真内容：

进行两组试验：第一组试验中，每次蒙特卡洛实验，采用本发明进行波达方向估计；第二组试验中，每次蒙特卡洛实验，采用L1-SVD算法进行波达方向检测，参照图2，为本发明和L1-SVD算法在波达方向均方根误差上的对比示意图。从图2可以看出，本发明的波达方向测量精度要高于L1-SVD算法，也就说明本发明在较低信噪比就可以达到较高的角度分辨能力。例如，当信噪比为5dB时，L1-SVD算法的波达方向均方根误差为7.5mrad，本发明的波达方向均方根误差仅为2mrad。参照图3，为利用本发明检测波达方向时的多径衰减系数估计误差随信噪比变化示意图。从图3可以看出，随着信噪比的增加，多径衰减系数的估计值越来越精确，说明：角度估计误差越小，定位精度越高，表明本发明是实际可行的。

综上所述，本发明可实现多径环境下高分辨的波达方向的估计，且方向估计精度高。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (5)

1.一种雷达用低空目标波达方向估计方法，其特征在于，包括以下步骤：

S1：雷达的天线阵列接收回波信号；对所述回波信号进行采样，得到数据矢量矩阵X；对数据矢量矩阵X进行奇异值分解，根据分解结果，生成对应的降维数据矢量X_c、以及噪声子空间矩阵E_n；

S2：根据雷达搜索角度范围的下限、雷达搜索角度范围的上限、以及雷达测角精度，得出角度采样点个数M、以及对应的第1个采样角度至第M个采样角度，计算与第i个采样角度θ_1i对应的第i镜像采样角度θ_2i，i取1至M；根据第i目标采样角度θ_1i和第i镜像采样角度θ_2i，得出对应的导向矢量阵A(θ_1i,θ_2i)；

S3：根据噪声子空间矩阵E_n、以及导向矢量阵A(θ_1i,θ_2i)，得出对应的多径衰减矢量构造复合导向矢量b(θ_1i)，构造感知矩阵B，B＝[b(θ₁₁),...,b(θ_1M)]；

S4：得出满足以下关系式的正则化参数μ：μ²≥Mσ²，σ²为所述回波信

号中的噪声信号的平均功率；建立如下范数约束式：

$\underset{{\mathrm{\Ψ}}_c}{min}{||{\mathrm{\Ψ}}_c||}_1$

$\begin{array}{cc}s.t.& {||X_c-{\mathrm{B\Ψ}}_c||}_2^2\≤{\mathrm{\μ}}^2\end{array}$

其中，Ψ_c为稀疏矢量，||·||₁表示L1范数，||·||₂表示L2范数，为：||Ψ_c||₁最小时的Ψ_c；然后，针对所述范数约束式，求解出Ψ_c；根据以下公式计算波达方向 其中，Ψ_c(i)为求解出的Ψ_c中的第i个元素，abs(Ψ_c(i))表示Ψ_c(i)的模值，表示min(abs(Ψ_c(i)))对应的θ_1i。

2.如权利要求1所述的一种雷达用低空目标波达方向估计方法，其特征在于，在步骤S1中，所述雷达的天线阵列包括N个等距离分布的阵元，每个阵元在接收到回波信号之后，进行L次采样；第n个阵元在进行第l次采样后得到的数据为x_n(l)，其中，n取1至N，l取1至L；然后构造数据矢量矩阵X，数据矢量矩阵X的行数为N，列数为L，数据矢量矩阵X的第n行第l列的元素为x_n(l)；

在步骤S1中，按照以下公式对数据矢量矩阵X进行奇异值分解：X＝UΛV^H，其中，U为N×N阶酉矩阵，Λ为N×L阶矩阵，V为L×L阶酉矩阵，V^H为V的共轭转置；然后根据以下公式生成对应的降维数据矢量X_c：

X_c＝UΛY＝XVY

其中，Y为N×1维向量，Y中第一行第一列元素为1，其余元素均为0；在N阶酉矩阵U中，取出后N-1列元素构成噪声子空间矩阵E_n。

3.如权利要求1所述的一种雷达用低空目标波达方向估计方法，其特征在于，在步骤S2中，根据以下公式计算得出角度采样点个数M：其中，α为雷达搜索角度范围的下限，β为雷达搜索角度范围的上限，Δθ为雷达测角精度；然后根据以下公式得出第i目标采样角度θ_1i：

θ_1i＝α+(i-1)·Δθ；

然后，根据以下公式计算与第i镜像采样角度θ_2i：

${\mathrm{\θ}}_{2i}=-arctan\frac{R{\mathrm{sin\θ}}_{1i}+2(h_a-h_s)}{R{\mathrm{cos\θ}}_{1i}},$

其中，R为目标与雷达之间的距离，h_a为雷达高度，h_s为多径效应的反射面高度，-h_s等于0；然后，构造目标导向列向量a(θ_1i)和镜像导向列向量a(θ_2i)，

$a\left({\mathrm{\θ}}_{1i}\right)=\left[\begin{array}{c}1\\ e^{-j\frac{2\mathrm{\π}f}{c}d{\mathrm{sin\θ}}_{1i}}\\ \·\\ \·\\ \·\\ e^{\text{-}j\frac{2\mathrm{\π}f}{c}(N-1)d{\mathrm{sin\θ}}_{1i}}\end{array}\right],a\left({\mathrm{\θ}}_{2i}\right)=\left[\begin{array}{c}1\\ e^{-j\frac{2\mathrm{\π}f}{c}d{\mathrm{sin\θ}}_{2i}}\\ \·\\ \·\\ \·\\ e^{\text{-}j\frac{2\mathrm{\π}f}{c}(N-1)d{\mathrm{sin\θ}}_{2i}}\end{array}\right]$

其中，f为雷达频率，d为雷达阵元间距，c为光速，则导向矢量阵A(θ_1i,θ_2i)为：A(θ_1i,θ_2i)＝[a(θ_1i)a(θ_2i)]。

4.如权利要求1所述的一种雷达用低空目标波达方向估计方法，其特征在于，在步骤S3中，按照以下公式计算噪声投影矩阵P： 表示E_n的共轭转置；按照以下公式得出多径衰减矢量

$\hat{p}({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})=\frac{{\[A^H({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})PA({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})\]}^{-1}w}{w^H{\[A^H({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})PA({\mathrm{\θ}}_{1i},{\mathrm{\θ}}_{2i})\]}^{-1}w}$

其中，w＝[1,0]^T。

5.如权利要求1所述的一种雷达用低空目标波达方向估计方法，其特征在于，在步骤S4中，针对所述范数约束式，利用凸优化设计工具箱CVX，求解出Ψ_c。