CN103018730B - 分布式子阵波达方向估计方法 - Google Patents

Abstract

本发明公开了一种基于分布式子阵的参差基线测角方法，主要解决分布式子阵对目标波达方向估计过程中出现角度模糊的问题，其实现过程是：从雷达回波中提取目标信号；使用数字波束形成法对目标信号进行角度粗测；对子阵进行划分，将每个子阵等分为两个新的子阵，形成四个新的子阵，产生四个基线，取其中的两个基线构造参差关系；对四个子阵进行波束合成，计算两个基线的相位差观测值；利用构造的参差关系计算模糊数的取值范围；引入代价函数，利用最小均方误差准则求得模糊数，并利用模糊数求得代价函数值，得到目标波达角。本发明采用构造两基线参差的处理方法来降低运算量，提高了对目标到达角估计的计算速度，可用于目标识别。

Description

分布式子阵波达方向估计方法

技术领域

本发明属于雷达信号处理技术领域，涉及分布式雷达测角，具体地说是针对分布式数字阵列雷达，提出一种基于分布式子阵的多分辨测角方法，可用于目标定位。

背景技术

波达方向DOA估计是阵列信号处理中的一个重要研究方向，在雷达、通信、声纳和导航等领域都得到了广泛的应用和研究。阵列的DOA估计性能受到阵列孔径的限制，阵列孔径越大，DOA估计性能越优，因此扩展孔径是提高测角精度的重要方法。早期微波合成的干涉雷达可以有效扩展孔径提高角度估计精度，但是过小的自由度限制了测角模糊问题的解决。为此Nilsson首次提出了分布式阵列雷达的概念，大大地扩展了阵列的自由度。分布式阵列在保持阵元数不变的基础上，通过扩大阵元间距来扩展阵列孔径，通过多个分布式小孔径阵列来提高测角精度与分辨率。但同时由于基线长度大于入射信号的半波长，出现相位缠绕，导致测角模糊。

由空域采样定理可知，对分布式阵列进行波束形成，将得到窄的主瓣及具有大量栅瓣的波束方向图，窄的主瓣提高了测角精度，而栅瓣导致了测角模糊。采用分布式孔径相参合成工作方式时，两个站的距离即基线长度可以适当的拉开；基线越长，合成阵列的等效孔径就越大，测角精度越高，但同时信号的相干性就越差、角度模糊越严重，所以需要研究在此种情况下两站或多站间协同探测方法、测角方法、角度解模糊方法，以满足雷达的需要。

为了得到精确的角度估计值，解决分布式阵列带来的角度模糊问题，目前主要采用以下几种方法：

（1）基于相位干涉仪的DOA估计：相位干涉仪解模糊主要通过不同的基线配置来实现，已有的方法包括：长短基线法，基于参差基线的相位差变化值法，阵列多组解模糊法，二次相位差解模糊法等。相位干涉仪具有测向精度高，结构简单，观测频带宽等优点，但是它的阵元利用率比较低，使得信噪比门限提高，而且对阵元的位置分布要求严格，需要满足特定条件，很大程度上限制了其在分布式阵列中的应用。

（2）基于双尺度ESPRIT的DOA估计：Zoltowski在《Direction finding with sparserectangular dual-size spatial invariance Array》中给出了双尺度ESPRIT解模糊算法。它是以分布式阵列的子阵获得精度低但无模糊的粗估计，使用子阵间的相位中心偏移来获得目标方向余弦的精估计，针对干涉阵带来的测角模糊问题，以粗估计为参考解精估计的模糊，从而得到高精度无模糊的角度估计。该方法和前面提到的长短基线法类似，但是需要进行特征值分解，解模糊时还需要进行方向余弦的配对，所以运算量极大，在快拍数较少的情况下，性能不是很好，而且随快拍数的增加算法性能增加不明显。

（3）基于子空间拟合的DOA估计：利用ESPRIT算法等对两个分布式子阵间的相位差进行估计，得到多个模糊的波达方向精估计值，然后再利用信号子空间拟合或者噪声子空间拟合关系，将这些可能值代入合适的代价函数中，选取其中最大或最小的估计值即为目标的DOA。该方法提供了比双尺度ESPRIT方法更高的角度估计精度，但是它不仅需要进行特征值分解，而且需要对MUSIC谱进行一维角度搜索，运算量较大，而且需要较多的快拍数，快拍数少时无法进行子空间拟合。

上述方法中，方法2和方法3能够得到高精度的DOA估计，但是需要特征值分解以及角度搜索，计算量大。方法1虽然避免了特征值分解，但是它的阵元利用率比较低，测角精度并不能达到系统设计的要求。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种分布式子阵波达方向估计方法，以进一步降低运算量，提高低信噪比情况下DOA的测角精度。

为实现上述目的，本发明的技术思路是:通过对分布式均匀子阵进行子阵划分构造出两基线的参差关系，计算产生的相位差，再利用余数定理进行解模糊处理，求得目标的DOA。具体实现步骤包括如下：

1）从雷达回波中提取目标信号：

X(t)＝A(μ)s(t)+n(t)；

其中，A(μ)为阵列对信号的导向矢量矩阵，s(t)为雷达发射信号，n(t)为均值为零、方差为1的高斯白噪声；

2）使用数字波束形成法DBF对目标信号X(t)进行角度粗测，得到目标信号的粗测角度θ_i0；

3）对分布式阵列中的两个子阵进行划分，将每个子阵等分为两个新的子阵，形成四个新的子阵，产生四个基线，取其中的两个基线D₁₃，D₁₄，分别产生两个相位差根据相位干涉仪原理，得到这两个相位差与两个基线长度的关系:

4）根据目标回波信号X(t)及粗测角度θ_i0，分别对四个子阵进行波束合成，将得到的合成后的信号平均求复角，计算出所取两个基线相位差的观测值和

5）根据基线长度大于半波长时会出现角度模糊的特性，将两个基线的模糊数分别设为k₁和k₂，利用余数定理计算模糊数k₁，k₂的取值范围：

$\left|k_1\right|\≤\min \{\frac{n_{1}d}{\mathrm{\λ}}+\frac{1}{2},p_1\}$

$\left|k_2\right|\≤\min \{\frac{n_{2}d}{\mathrm{\λ}}+\frac{1}{2},p_2\},$

其中，d为阵元间距，λ表示波长，n₁，n₂分别为两个基线长度与阵元间距的比值，n₁＝F·p₁，n₂＝F·p₂，F为n₁，n₂的最大公约数，p₁，p₂为互素的正整数；

6）引入代价函数L，令：

L₁＝n₂k₁+h₁

L₂＝n₁k₂+h₂,

其中，L₁，L₂分别为两个基线的代价函数，h₁，h₂为考虑噪声影响时L₁，L₂的测向误差；

7）利用最小均方误差准则，在上述取值范围|k₁|，|k₂|内搜索，得到使得(L₁-L₂)²最小的一组模糊数k₁,k₂，即为所求模糊数；

8）将上述两个模糊数k₁,k₂分别代入代价函数L₁，L₂的表达式，求得代价函数L₁，L₂的值，得到目标波达角： $\mathrm{\θ}=\arcsin \left(\frac{(L_1+L_2)\mathrm{\λ}}{{2n}_1{n}_{2}d}\right).$

本发明与现有技术相比具有如下优点：

（1）现有的分布式子阵的波达方向估计方法在解模糊时多采用MUSIC、ESPRIT等超分辨算法，需要进行特征值分解以及角度的多维搜索，运算量比较大，而本发明采用构造两基线参差的处理方法，不需要进行角度搜索，运算量大大降低，进而提高了运算速度。

（2）现有的分布式子阵的波达方向估计方法在快拍数较低的情况下，进行角度搜索时不够精确，而本发明采用DBF粗测的角度进行子阵阵元的波束形成，降低了角度解模糊的SNR门限，提高了测角精度。

附图说明

图1是本发明的流程图；

图2是本发明中分布式子阵的旋转不变子阵划分的示意图；

图3是本发明中正确解模糊概率和信噪比的关系示意图；

图4是用本发明与基于双尺度ESPRIT的DOA估计算法对目标DOA估计的性能比较图。

具体实施方式

下面结合附图详细说明本发明的内容和效果。

参照图1，本发明包括如下步骤：

步骤1：从雷达回波中提取目标信号。

本发明中分布式子阵的旋转不变子阵划分模型如图2所示。图2中分布式阵列由两个基线距离为Δ_s的子阵组成，每个子阵内等间隔分布N_s个各向同性阵元，因此总阵元数为N＝2N_s，阵元间距为d≤λ/2，λ为信号波长，子阵间距远大于阵元间距。

假设窄带远场信号以角度θ入射到阵列上，以第一个子阵的中心阵元为参考点，则参考子阵对信号的导向矢量为：

$a_{\mathrm{sub}}\left(\mathrm{\μ}\right)={[e^{-j(N_s-1)/{2\mathrm{\ω}}_k},e^{-j(N_s-3)/{2\mathrm{\ω}}_k},...,e^{j(N_s-1)/{2\mathrm{\ω}}_k}]}^T,$

其中ω_k＝(2π/λ)·d·μ是阵元间距d对信号的相移，μ＝sinθ是方向余弦，上标T表示转置，因此整个阵列对信号的导向矢量为：

$A\left(\mathrm{\μ}\right)=a_{\mathrm{sd}}\left(\mathrm{\μ}\right)\⊗a_{\mathrm{sub}}\left(\mathrm{\μ}\right),$

其中a_sd(μ)表示Δ_s引起的子阵间的相移向量，

$a_{\mathrm{sd}}\left(\mathrm{\μ}\right)={[e^{-j(N_s-1)/{{2\mathrm{\ω}}_k}^f},e^{-j(N_s-3)/{{2\mathrm{\ω}}_k}^f},...,e^{j(N_s-1)/{{2\mathrm{\ω}}_k}^f}]}^T,$

其中ωk_k ^f＝(2π/λ)·Δ_s·μ是子阵间距Δ_s对信号的相移。

从图2信号模型中得到阵元接收的目标信号：

X(t)＝A(μ)s(t)+n(t)，

其中，s(t)为雷达发射信号，n(t)为均值为零、方差为1的高斯白噪声。

步骤2：使用数字波束形成法DBF对目标信号进行角度粗测，得到目标信号的粗测角度θ_i0：

${\mathrm{\θ}}_{i0}=\arg \underset{\mathrm{\θ}}{\max }\left(\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\left|a^H\left(\mathrm{\θ}\right)X\left(t_l\right)\right|}^2\right),$

其中：argmax为寻找具有最大代价函数的参量，L表示快拍数，a(θ)＝[e^{jκ°·sin(θ)},e^{jκ·1·sin(θ)},…,e^{jκ(N-1)sin(θ)}]^T，θ表示目标搜索角度，κ表示波数，N表示阵元个数，X(t_l)表示阵元接收到的目标信号，上标T表示转置，上标H表示共轭转置。

步骤3：对分布式阵列中的两个子阵进行划分，得到两个相位差与两个基线长度的关系。

本发明中分布式子阵的旋转不变子阵的划分如图2所示。

3a）将每个子阵等分为两个新的子阵，即图中所示的X₁~X₄，形成四个新的子阵；

3b）根据这四个子阵都具有旋转不变关系，写出对应的四组基线长度分别为：

D₁₂＝D₃₄＝(N_s-N_s/2)d＝N_s/2·d

D₁₃＝D₂₄＝(N_d+N_s)d

D₁₄＝(N_s+N_d+N_s-N_s/2)d＝(N_d+3N_s/2)d

D₂₃＝(N_d+N_s/2)d，

其中，N_d为阵列间隔阵元数，D₁₂为第一个子阵与第二个子阵的基线长度，D₁₃第一个子阵与第三个子阵的基线长度，D₁₄第一个子阵与第四个子阵的基线长度，D₂₃为第二个子阵与第三个子阵的基线长度，D₂₄为第二个子阵与第四个子阵的基线长度，D₃₄为第三个子阵与第四个子阵的基线长度；

3c）取其中的两个基线D₁₃，D₁₄，根据相位干涉仪原理分别产生两个相位差

3d）得出两个相位差与两个基线长度的关系为：

步骤4：根据目标回波信号X(t)及粗测角度θ_i0，分别对四个子阵进行波束合成，将得到的合成后的信号平均求复角，计算出所取两个基线相位差的观测值和

根据目标回波信号X(t)及粗测角度θ_i0，分别对四个子阵进行波束合成，得到合成后的信号分别为X₁(t)，X₂(t)，X₃(t)，X₄(t)，将（X₁(t)，X₃(t)）和（X₂(t)，X₄(t)）平均求复角，可以得到相位差的观测值同理，对（X₁(t)，X₄(t)）求复角可以得到相位差的观测值

步骤5：根据基线长度大于半波长时会出现角度模糊的特性，将两个基线的模糊数分别设为k₁和k₂，利用余数定理计算模糊数k₁，k₂的取值范围。

5a）根据基线长度大于半波长会出现角度模糊和相位差是以2π为周期模糊的特性，得到相位差与其观测值的以下关系：

5b）设n₁，n₂分别为两个基线长度与阵元间距的比值，则两个基线长度可以表示为D₁₃＝n₁·d，D₁₄＝n₂·d；根据步骤3得到的相位差和基线长度的关系，进而可以得到相位差和n₁，n₂的关系：

5c）联立步骤5a，5b中的三个方程，在理想无噪声的情况下，得到关系式：

上式是一个除数为整数的实数域内的同余方程，k₁，k₂在一定范围内存在唯一解。

5d）根据n₁，n₂分别是两个基线长度与阵元间距比值的定义，得到n₁＝N_d+N_s＝F·p₁，n₂＝N_d+3N_s/2＝F·p₂，其中，F为n₁，n₂的最大公约数，p₁，p₂为互素的正整数；

对于步骤5c中的同余方程，由余数定理可以得到模糊数k₁，k₂取值范围的一个约束条件：

-p₁≤k₁≤p₁

-p₂≤k₂≤p₂；

5e）给出模糊数k₁，k₂取值范围的另一个约束条件，即由前面的关系式D₁₃＝n₁·d，将k₁表示为：

根据sinθ最大为1，得到k₁的范围：

$\left|k_1\right|\≤\frac{n_{1}d}{\mathrm{\λ}}\·1-\frac{-\mathrm{\π}}{2\mathrm{\π}}=\frac{n_{1}d}{\mathrm{\λ}}+\frac{1}{2},$

同理可得k₂的范围：

$\left|k_2\right|\≤\frac{n_{2}d}{\mathrm{\λ}}\·1-\frac{-\mathrm{\π}}{2\mathrm{\π}}=\frac{n_{2}d}{\mathrm{\λ}}+\frac{1}{2};$

5f）综合上面给出的k₁，k₂取值范围的两个约束条件，得到k₁，k₂的取值范围为：

$\left|k_1\right|\≤\min \{\frac{n_{1}d}{\mathrm{\λ}}+\frac{1}{2},p_1\}$

$\left|k_2\right|\≤\min \{\frac{n_{2}d}{\mathrm{\λ}}+\frac{1}{2},p_2\}.$

步骤6：引入代价函数L。

当考虑噪声的影响时，相位差存在误差，上述的关系式不再成立，为了求得模糊数，引入代价函数L，令：

L₁＝n₂k₁+h₁

L₂＝n₁k₂+h₂,

其中，L₁，L₂分别为两个基线的代价函数，h₁，h₂为考虑噪声影响时L₁，L₂的测向误差， 为噪声引起的相位误差。

步骤7：利用最小均方误差准则，在上述取值范围|k₁|，|k₂|内搜索，得到使得(L₁-L₂)²最小的一组模糊数k₁,k₂，即为所求模糊数。

步骤8：利用这两个模糊数k₁,k₂求得目标波达角。

将k₁，k₂分别代入代价函数L₁，L₂的表达式，即可求得代价函数L₁，L₂，根据目标波达角θ与代价函数L₁，L₂的关系：最后得到目标波达角：

$\mathrm{\θ}=\arcsin \left(\frac{(L_1+L_2)\mathrm{\λ}}{{2n}_1{n}_{2}d}\right).$

本发明的效果通过以下计算仿真进一步说明：

1.仿真条件

仿真条件为以下雷达参数：子阵阵元数N_s＝16，阵列间隔阵元数N_d＝80，阵元间隔为半波长，目标个数P＝1，快拍数L＝1，角度搜索个数n＝180。

2.仿真内容

仿真1，用本发明对分布式子阵的正确解模糊概率进行仿真，仿真结果如图3所示。其中横轴表示目标信噪比从0dB至30dB变化，纵轴表示正确解模糊概率。从图3中可以看出，当信噪比大于8dB时，本发明可以完全正确解模糊。

仿真2，用现有的基于双尺度ESPRIT的DOA估计算法和本发明分别对目标进行测角精度仿真，仿真选取的目标参数：目标仰角10度，蒙特卡罗实验次数100次，仿真结果如图4所示，其中横轴表示信噪比从0dB至30dB变化，纵轴表示测角误差。图4中双尺度ESPRIT表示基于双尺度ESPRIT的DOA估计算法在信噪比按照横轴变化时的测角均方误差，参差基线表示本发明在信噪比按照横轴变化时的测角均方误差。从图4可以得出，在低信噪比的情况下，本发明要比基于双尺度ESPRIT的DOA估计算法的测角误差小。

仿真3，计算复杂度比较。表1列出了本发明和现有双尺度ESPRIT算法的计算复杂度。

表1本发明和现有算法的计算复杂度

从表1中可以看出，本发明的计算复杂度为O{3.03×10³}，而双尺度ESPRIT算法的计算复杂度为O{3.38×10⁴}，本发明的计算复杂度降低了一个数量级。而且Ns越大，阵元数越多，算法复杂度降低越多。可见本发明能大大降低雷达信号处理的运算量，提高目标波达角估计的计算速度。

Claims (3)

1.一种分布式子阵波达方向估计方法，其特征在于包括以下步骤：

1)从雷达回波中提取目标信号:

X(t)＝A(μ)s(t)+n(t)；

其中，A(μ)为阵列对信号的导向矢量矩阵，s(t)为雷达发射信号，n(t)为均值为零、方差为1的高斯白噪声；

2)使用数字波束形成法DBF对目标信号X(t)进行角度粗测，得到目标信号的粗测角度θ_i0；

3)对分布式阵列中的两个子阵进行划分，将每个子阵等分为两个新的子阵，且两个新的子阵具有重叠阵元，即将第一子阵化分为两个新子阵X1和X2，将第二子阵化分为两个新子阵X3和X4，形成四个新的子阵，产生四个基线，取其中的两个基线，分别产生两个相位差根据相位干涉仪原理，得到这两个相位差与两个基线的长度关系:

其中，D₁₃为第一个新的子阵X1与第三个新的子阵X3的基线长度，D₁₄为第一个新的子阵X1与第四个新的子阵X4的基线长度；

4)根据目标回波信号X(t)及粗测角度θ_i0，分别对四个新的子阵进行波束合成，将得到的合成后的信号平均求复角，计算出所取两个基线相位差的观测值和

5)根据基线长度大于半波长时会出现角度模糊的特性，将两个基线的模糊数分别设为k₁和k₂，利用余数定理计算模糊数k₁，k₂的取值范围:

$\left|k_1\right|\≤\min \{\frac{n_{1}d}{\mathrm{\λ}}+\frac{1}{2},p_1\}$

$\left|k_2\right|\≤\min \{\frac{n_{2}d}{\mathrm{\λ}}+\frac{1}{2},p_2\},$

其中，d为阵元间距，λ表示信号波长，d≤λ/2，n₁，n₂分别为两个基线长度与阵元间距的比值，n₁＝F·p₁，n₂＝F·p₂，F为n₁，n₂的最大公约数，p₁，p₂为互素的正整数；

6)引入代价函数L，令：

L₁＝n₂k₁+h₁

L₂＝n₁k₂+h₂,

其中，L₁，L₂分别为两个基线的代价函数，h₁，h₂为考虑噪声影响时L₁，L₂的测向误差；

7)利用最小均方误差准则，在上述取值范围|k₁|，|k₂|内搜索，得到使得(L₁-L₂)²最小的一组模糊数k₁,k₂，即为所求模糊数；

8)将上述两个模糊数k₁,k₂分别代入代价函数L₁，L₂的表达式，求得代价函数L₁，L₂的值，得到目标波达角:

2.根据权利要求1所述的分布式子阵波达方向估计方法，其中步骤2)所述的使用数字波束形成法DBF对目标信号进行角度粗测，是通过如下公式进行:

${\mathrm{\θ}}_{i0}=\arg \underset{\mathrm{\θ}}{\max }\left(\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\left|a^H\left(\mathrm{\θ}\right)X\left(t_l\right)\right|}^2\right)$

其中：θ_i0为目标信号的粗测角度，arg max为寻找具有最大代价函数的参量，L表示快拍数，a(θ)＝[e^{jκ·0·sin(θ)},e^{jκ·1·sin(θ)},…,e^{jκ·(N-1)·sin(θ)}]^T，θ表示目标搜索角度，κ表示波数，N表示阵元个数，X(t_l)表示阵元接收到的目标信号，上标T表示转置，上标H表示共轭转置。

3.根据权利要求1所述的分布式子阵波达方向估计方法，其中步骤3)所述的对分布式阵列中的两个子阵进行划分，按如下步骤进行：

3a)将每个子阵等分为两个新的子阵，形成四个新的子阵；

3b)根据这四个新的子阵都具有旋转不变关系，写出对应的四组基线长度分别为:

D₁₂＝D₃₄＝(N_s-N_s/2)d＝N_s/2·d

D₁₃＝D₂₄＝(N_d+N_s)d

D₁₄＝(N_s+N_d+N_s-N_s/2)d＝(N_d+3N_s/2)d

D₂₃＝(N_d+N_s/2)d，

其中，N_s为单个子阵所包含的阵元数，N_d为阵列间隔阵元数，d为阵元间距，D₁₂为第一个新的子阵X1与第二个新的子阵X2的基线长度，D₁₃为第一个新的子阵X1与第三个新的子阵X3的基线长度，D₁₄为第一个新的子阵X1与第四个新的子阵X4的基线长度，D₂₃为第二个新的子阵X2与第三个新的子阵X3的基线长度，D₂₄为第二个新的子阵X2与第四个新的子阵X4的基线长度，D₃₄为第三个新的子阵X3与第四个新的子阵X4的基线长度。