CN103399308A - 主瓣和旁瓣干扰背景下雷达目标角度快速估计方法 - Google Patents

Abstract

本发明公开了一种主瓣和旁瓣干扰背景下雷达目标角度快速估计方法，由雷达回波数据得到含有目标角度信息的复数据列矢量；计算把扩展导向矢量变换为多项式的变换矩阵；把雷达系统归一化空间频率的似然函数的导数转化为4组低阶多项式；求4组低阶多项式的实根；在实际波束宽度边缘点和多项式实根中求出目标归一化空间频率，得到目标角度。本发明使归一化空间频率的似然函数能够直接求导，其中似然函数的极点由低阶实多项式的根确定，克服了现有技术计算复杂度高、数值稳定性差和不适用于大阵列的不足，具有计算复杂度低，数值稳定性好的优点，特别适用于大阵列在主瓣和旁瓣干扰背景下雷达目标角度的快速估计。

Description

主瓣和旁瓣干扰背景下雷达目标角度快速估计方法

技术领域

本发明属于雷达测角技术领域，主要涉及测角中的信号处理，具体是一种主瓣和旁瓣干扰背景下雷达目标角度快速估计方法，可用于自适应多通道雷达在干扰环境中的测角。

背景技术

测量目标的角度是现代侦察跟踪雷达的重要功能。1976年，Davis等人采用网格搜索最大似然方法得到自适应最优角度估计，这种方法的优点是能自适应抑制干扰，但是这种方法需要在波束宽度内取非常多的角度网格搜索出似然函数的最大值，测角准确度和用于搜索的角度网格数有关，网格数越多，准确度越高，所以这种方法的实现需要非常高的计算复杂度。

2010年，吴建新等人提出了一种传统求根最大似然方法进行雷达目标测角。这种方法把主波束内的导向矢量矩阵用保角变换，转化为多项式系数矩阵，从而把似然函数的导数转化为一个4N阶实多项式，其中N是阵元数，似然函数的最大值可以通过在4N阶多项式的实根和波束宽度的边缘点中求得，似然函数最大值对应的角度值就为雷达目标的角度。如果阵元N很小，这种方法的计算优势是很明显的，并且可以保证数值稳定性。但是，随着阵元N的增加，由于高阶多项式求根计算复杂而且数值稳定性较差，导致雷达测角的过程时间长、效率低，这些优势就减小了。以这种方法对于大阵列雷达来说甚至无法实现满足实战需要的目标测角。

发明内容

本发明的目的在于克服现有技术的不足，提出一种计算复杂度低、数值稳定性好、与阵元数无关的主瓣和旁瓣干扰背景下雷达目标角度快速估计方法。

本发明是一种主瓣和旁瓣干扰背景下雷达目标角度快速估计方法，实现本发明的基本思路是：首先利用虚拟波束内扩展导向矢量组成的扩展导向矢量矩阵的低秩特性，把虚拟波束内的任意一个扩展导向矢量分解为变换矩阵和多项式矢量乘积的形式，从而对归一化空间频率的似然函数求导用4组低阶实多项式表示，并且似然函数的极点由4组低阶实多项式的根确定。最后，在实际波束宽度边缘点和4组多项式的实根中求出似然函数的最大值，最大值对应的归一化空间频率就为目标归一化空间频率的最大似然估计。

本发明雷达目标角度快速估计过程包括如下步骤：

步骤1：雷达通过天线接收回波数据，该数据按阵元通道和距离单元二维矩阵的格式存储，回波数据中包括主瓣干扰和旁瓣干扰、目标信号和接收机噪声；

步骤2：用雷达回波数据估计干扰和噪声的协方差矩阵，用自适应滤波、单元平均恒虚警检测方法确定目标所在的距离单元，用干扰和噪声的协方差矩阵的逆矩阵和目标所在距离单元的雷达回波数据构造白化数据矩阵A，用该白化数据矩阵信息和所述协方差矩阵的逆矩阵信息构造两个列矢量，利用两个列矢量信息得到含有目标角度信息的复数据列矢量c；

步骤3：计算把扩展导向矢量变换为多项式的变换矩阵F；

(3a)求雷达虚拟波束宽度内扩展导向矢量的低秩子空间U_P；

(3b)在虚拟波束宽度内，把不同角度的扩展导向矢量在子空间的投影系数矢量用多项式进行拟合，得到系数矩阵B；

(3c)用低秩子空间U_P、系数矩阵B和对角阵D相乘，得到把扩展导向矢量转化为多项式的变换矩阵：

F＝U_PDB

其中D是大小为P×P的对角阵，P是低秩子空间U_P的维数，D＝diag([ε₁ ε₂ … ε_P])，diag表示对角阵，D的对角线元素ε_p满足

i为虚数单位；

步骤4：把实际波束宽度等分成4部分，每一部分的波束宽度小于虚拟波束宽度，运用含有目标角度信息的复数据列矢量c、变换矩阵F和与实际波束宽度等分有关的对角阵Λ_m，把雷达系统归一化空间频率v的似然函数J(v)的导数J(v)′转化为4组低阶多项式J_m(v)′＝β(v)^HF^HΛ_mc，m＝1，2，3，4，

其中上标′表示求导，β(v)＝[1 v … v^Q]^T为多项式矢量，Q为多项式阶数，上标T表示转置，上标H表示共轭转置，

diag表示对角阵， $s_1\left(v\right)=e^{-i2\mathrm{\π}(2N-2)v}{\left[\begin{array}{cccc}1& e^{i2\mathrm{\πv}}& ...& e^{i2\mathrm{\π}(4N-4)v}\end{array}\right]}^T$ 是大小为4N-3的扩展导向矢量，i是虚数单位，v₀是视线方向归一化空间频率，

是实际波束宽度，N是阵元数，归一化空间频率v的取值范围是 $v\∈[v_0-\frac{1}{2}{\mathrm{\Δ}}_N,v_0+\frac{1}{2}{\mathrm{\Δ}}_N];$

步骤5：用牛顿法求4组低阶多项式的实根；

步骤6：把实际波束宽度的边缘点和4组低阶多项式的实根代入归一化空间频率似然函数J(v)，计算出似然函数的最大值，最大值对应的归一化空间频率就为目标的归一化空间频率的最大似然估计，用归一化空间频率和角度的转换关系，求出目标的角度。

本发明针对现有技术中网格搜索最大似然方法计算复杂度高和传统的求根最大似然方法对于大阵列计算复杂度高、数值稳定性差等诸多不足，本发明通过把雷达虚拟波束宽度内扩展导向矢量用变换矩阵转化为多项式，从而使归一化空间频率的似然函数的极值从4组低阶多项式的根中求得，快速实现了自适应阵列最大似然角度估计，既降低了计算复杂度，又确保了数值稳定性，更重要的是多项式的阶数和阵元数是独立的，因此对大阵列是非常有益的，测角效率高，性能稳定。

本发明的实现还在于：步骤2中构造含有目标角度信息的复数据列矢量c的方法包括：

(2a)利用干扰和噪声的协方差矩阵的逆矩阵R^-1和目标所在距离单元的雷达回波数x_t据构造白化数据矩阵A，利用白化数据矩阵的对角元素构造列矢量1；

白化数据矩阵A按下式构造：

$A=R^{-1}{x}_t{x}_t^H{R}^{-1}$

其中上标-1表示矩阵求逆，上标H表示共轭转置，

利用白化数据矩阵A的对角线元素，按下两式构造列矢量1：

$a_k=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}A(N-k+n,n),$ k＝1，2，...，N

$a_1={\left[\begin{array}{ccccccc}{a}_1& \·\·\·& a_{N-1}& a_N& a_{N-1}^{*}& \·\·\·& a_1^{*}\end{array}\right]}^T$

其中a₁为大小为(2N-1)×1的列矢量1，上标*表示共轭，上标T表示转置，a_k是构造a₁的基本元素，k的取值范围为1，2，…，N，N为阵元数；

(2b)利用干扰和噪声的协方差矩阵的逆矩阵R^-1，按下两式构造列矢量2：

$r_k=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{R}^{-1}(N-k+n,n),$ k＝1，2，...，N

$a_2={\left[\begin{array}{ccccccc}{r}_1& \·\·\·& r_{N-1}& r_N& r_{N-1}^{*}& \·\·\·& r_1^{*}\end{array}\right]}^T$

其中a₂为大小为(2N-1)×1的列矢量2，r_k是构造a₂的基本元素，k的取值范围为1，2，…，N；

(2c)用列矢量1、列矢量2和对角阵D₁按下式求得一个新矩阵：

$C=D_1^H{a}_1{a}_2^H+a_1{a}_2^H{D}_1$

其中C是大小为(2N-1)×(2N-1)的新矩阵，D₁是对角线元素是首项为-i2π(N-1)、公差是i2π的等差数列的对角阵，即D₁＝diag([-i2π(N-1)…0…i2π(N-1)])，diag表示对角阵，i为虚数单位。

(2d)利用新矩阵C，按下两式构造出含有目标角度信息的复数据列矢量c：

$c_k=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}C(N-k+n,n),$ k＝1，2，...，2N-1

$c={\left[\begin{array}{ccccccc}{c}_1& \·\·\·& c_{2N-2}& c_{2N-1}& c_{2N-2}^{*}& \·\·\·& c_1^{*}\end{array}\right]}^T$

其中c为大小为(4N-3)×1的含有目标角度信息的复数据列矢量，c_k是构造c的基本元素，k的取值范围为1，2，…，2N-1。

本发明针对现有技术中网格搜索最大似然方法不能对归一化空间频率似然函数直接求导，只能采用在波束宽度内进行网格搜索，计算复杂度高的不足，对回波数据进行构造含有目标角度信息的复数据列矢量的处理，使直接对归一化空间频率似然函数求导成为可能，因而不用一一搜索目标角度，只用在似然函数的导数的实根和实际波束宽度的边缘点处搜索，就能求出目标角度，大大降低了计算复杂度。

本发明的实现还在于：步骤(3a)中得到雷达虚拟波束宽度内扩展导向矢量的低秩子空间U_P的过程包括：

(3a1)雷达天线阵元数为N，虚拟波束宽度为

在虚拟波束宽度内构造sinc矩阵，记为G，它的每一个元素按下式计算：

$g_{\mathrm{xy}}\left(\mathrm{\Δ}\right)=\mathrm{\Δ}{e}^{i2\mathrm{\π}(x-y)v_0}\sin c\left(\mathrm{\Δ}(x-y)\right),$ x＝1，2，…，4N-3；y＝1，2，…，4N-3；

(3a2)对sinc矩阵G做特征分解，选取大的特征值对应的特征向量构成U_P，选取的特征向量的个数为6，即子空间的维数P为6。

本发明针对传统的求根最大似然方法中对于大阵列计算复杂度高、数值稳定性差的不足，把虚拟波束宽度内扩展导向矢量矩阵用低秩子空间表示，低秩子空间拟合误差对于不同的阵元数是基本相同的，表明张成虚拟波束空间的基矢量数量和阵元数是独立的，这对于减少计算复杂度来说是非常有用的性质。

本发明的实现还在于：步骤(3b)中把虚拟波束宽度内，不同角度的扩展导向矢量在低秩子空间的投影系数矢量用多项式进行拟合，得到系数矩阵B的过程包括有：

(3b1)在虚拟波束宽度内取不同角度的样本点v₁，v₂，...，v_J，其中J是选取角度样本点的个数；

(3b2)按下式构造不同角度对应的多项式矩阵v：

V＝[β(v₁) β(v₂) … β(v_J)]

其中 $\mathrm{\β}\left(v_j\right)={\left[\begin{array}{cccc}1& v_j& ...& v_j^Q\end{array}\right]}^T,$ j＝1，2，...，J是大小为(Q+1)×1的实矢量，Q是多项式的阶数，v的大小为(Q+1)×J；

(3b3)按下式计算不同角度的扩展导向矢量在子空间每一维矢量的投影系数a_p(v_j)，

p＝1，...，P，j＝1，...，J，

其中u_p是低秩子空间U_P的第p维列矢量， $s_1\left(v_j\right)=e^{-i2\mathrm{\π}(2N-2)v_j}{\left[\begin{array}{cccc}1& e^{i2\mathrm{\π}{v}_j}& ...& e^{i2\mathrm{\π}(4N-4)v_j}\end{array}\right]}^T,$ j＝1，2，...，J是大小为(4N-3)×1的扩展导向矢量；

把不同维的投影系数按不同角度排成一列，得到大小为J×1投影系数矢量a_p，a_p＝[a_p(v₁) a_p(v₂) … a_p(v_J)]^T，p＝1，...，P；

(3b4)按下式计算子空间每一维矢量u_p的投影系数矢量a_p对应的多项式拟合系数矢量b_p：

b_p＝(VV^T)^-1Va_p，p＝1，...，P，

其中b_p的大小是(Q+1)×1，上标-1表示矩阵求逆，

把所有维的多项式拟合系数矢量排成矩阵，并求转置，得到大小为P×(Q+1)系数矩阵B：

B＝[b₁ b₂ …b_P]^T。

本发明针对传统的求根最大似然角度估计方法中对于大阵列计算复杂度高、数值稳定性差的不足，把虚拟波束宽度内扩展导向矢量在低秩子空间的投影系数矢量，用多项式进行拟合，多项式阶数较低时，多项式拟合误差非常小，表明虚拟波束宽度内扩展导向矢量在低秩子空间的投影系数矢量可以用低阶多项式很好地近似，并且多项式拟合误差对于不同的阵元数是基本相同的，表明多项式的阶数和阵元数是独立的，这对于减少计算复杂度来说是非常有用的性质。

与现有技术相比，本发明的技术优势是：

第一，本发明通过对回波数据进行构造含有目标角度信息的复数据列矢量和利用虚拟波束宽度内扩展导向矢量组成的扩展导向矢量矩阵的低秩特性，把虚拟波束宽度内任意一个扩展导向矢量分解为变换矩阵和多项式矢量乘积形式的处理，使对归一化空间频率的似然函数直接求导成为可能，并且把归一化空间频率的似然函数的导数转化为多项式，降低了计算复杂度。

第二，本发明通过把实际波束宽度等分成4部分，将归一化空间频率的似然函数求导转化的多项式表示为4组低阶实多项式，从而归一化空间频率的似然函数的最大值可以从实际波束宽度边缘点和4组低阶实多项式的根中确定，最大值对应的归一化空间频率就为目标归一化空间频率的最大似然估计，由于多项式的阶数降低且与阵元数无关，因此这种技术具有计算复杂度低，数值稳定性好，适用于大阵列的优点。

附图说明

图1是本发明的流程图；

图2是本发明仿真的对于不同的阵元数，子空间拟合误差随着子空间维数变化的曲线图；

图3是本发明仿真的对于不同的阵元数，多项式拟合误差随着多项式阶数变化的曲线图；

图4是本发明仿真的在实际波束宽度内网格搜索最大似然方法似然函数值的曲线图及传统求根最大似然方法似然函数的极值点和本发明提出的求根最大似然方法似然函数的极值点图；

图5是本发明与网格搜索最大似然方法、传统求根最大似然方法针对目标角度和均方根误差的曲线图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细描述：

实施例1

本发明是一种主瓣和旁瓣干扰背景下雷达目标角度快速估计方法，地面雷达或机载雷达工作时，接收回波数据进行实时处理，估计测量出目标的角度，现有技术计算复杂度高，实时处理速度较慢，本发明针对这种现状提出一种对于大阵列雷达系统的快速测角新方法，该方法是在主瓣和旁瓣均存在干扰背景下实现的，无论是地面雷达或机载雷达，本发明都可应用，参见图1，本发明对雷达目标角度的快速估计包括以下步骤：

步骤1：雷达通过天线接收回波数据，雷达天线为N个阵元的均匀线性阵列，回波数据包括目标信号、干扰信号和接收机噪声，干扰信号包括主瓣干扰和旁瓣干扰，该数据按阵元通道和距离单元二维矩阵的格式存储。

步骤2：由雷达回波数据经过一系列处理得到含有目标角度信息的复数据列矢量c，即用雷达回波数据估计干扰和噪声的协方差矩阵R，用自适应滤波、单元平均恒虚警检测方法确定目标所在的距离单元x_t，用干扰和噪声的协方差矩阵的逆矩阵R^-1和目标所在距离单元的雷达回波数据x_t构造白化数据矩阵A，用该白化数据矩阵信息和所述协方差矩阵的逆矩阵信息构造两个列矢量，利用两个列矢量信息得到含有目标角度信息的复数据列矢量c。本发明通过对回波数据进行构造含有目标角度信息的复数据列矢量的处理，使对目标归一化空间频率的似然函数直接求导成为可能。

构造含有目标角度信息的复数据列矢量c的过程包括：

(2a)利用干扰和噪声的协方差矩阵的逆矩阵R^-1和目标所在距离单元的雷达回波数x_t据构造白化数据矩阵A，利用白化数据矩阵的对角元素构造列矢量1；

白化数据矩阵A按下式构造：

$A=R^{-1}{x}_t{x}_t^H{R}^{-1}$

其中上标-1表示矩阵求逆，上标H表示共轭转置，

利用白化数据矩阵A的对角线元素，按下两式构造列矢量1：

$a_k=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}A(N-k+n,n),$ k＝1，2，...，N

$a_1={\left[\begin{array}{ccccccc}{a}_1& \·\·\·& a_{N-1}& a_N& a_{N-1}^{*}& \·\·\·& a_1^{*}\end{array}\right]}^T$

其中a₁为大小为(2N-1)×1的列矢量1，上标*表示共轭，上标T表示转置，a_k是构造a₁的基本元素，k的取值范围为1，2，…，N，N为阵元数；

(2b)利用干扰和噪声的协方差矩阵的逆矩阵R^-1，按下两式构造列矢量2：

$r_k=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{R}^{-1}(N-k+n,n),$ k＝1，2，...，N

$a_2={\left[\begin{array}{ccccccc}{r}_1& \·\·\·& r_{N-1}& r_N& r_{N-1}^{*}& \·\·\·& r_1^{*}\end{array}\right]}^T$

其中a₂为大小为(2N-1)×1的列矢量2，r_k是构造a₂的基本元素，k的取值范围为1，2，…，N；

(2c)用列矢量1、列矢量2和对角阵D₁按下式求得一个新矩阵：

$C=D_1^H{a}_1{a}_2^H+a_1{a}_2^H{D}_1$

其中C是大小为(2N-1)×(2N-1)的新矩阵，D₁是对角线元素是首项为-i2π(N-1)、公差是i2π的等差数列的对角阵，即D₁＝diag([-i2π(N-1) … 0 … i2π(N-1)])，diag表示对角阵，i为虚数单位。

(2d)利用新矩阵C，按下两式构造出含有目标角度信息的复数据列矢量c：

$c_k=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}C(N-k+n,n),$ k＝1，2，...，2N-1

$c={\left[\begin{array}{ccccccc}{c}_1& \·\·\·& c_{2N-2}& c_{2N-1}& c_{2N-2}^{*}& \·\·\·& c_1^{*}\end{array}\right]}^T$

其中c为大小为(4N-3)×1的含有目标角度信息的复数据列矢量，c_k是构造c的基本元素，k的取值范围为1，2，…，2N-1。

步骤3：为了把虚拟波束宽度内，任意扩展导向矢量转化为多项式，计算变换矩阵F。

(3a)求雷达虚拟波束宽度内扩展导向矢量的低秩子空间U_P。

(3a1)雷达天线阵元数为N，虚拟波束宽度为

在虚拟波束宽度内构造sinc矩阵，记为G，它的每一个元素按下式计算：

$g_{\mathrm{xy}}\left(\mathrm{\Δ}\right)=\mathrm{\Δ}{e}^{i2\mathrm{\π}(x-y)v_0}\sin c\left(\mathrm{\Δ}(x-y)\right),$ x＝1，2，…，4N-3；y＝1，2，…，4N-3；

(3a2)对sinc矩阵G做特征分解，选取大的特征值对应的特征向量构成U_P，选取的特征向量的个数为6，即子空间的维数P为6。

(3b)在虚拟波束宽度内，把不同角度的扩展导向矢量在子空间的投影系数矢量用多项式进行拟合，得到系数矩阵B。

(3b1)在虚拟波束宽度内取不同角度的样本点v₁，v₂，...，v_J，其中J是选取角度样本点的个数；

(3b2)按下式构造不同角度对应的多项式矩阵v：

V＝[β(v₁ )β(v₂) … β(v_J)]

其中 $\mathrm{\β}\left(v_j\right)={\left[\begin{array}{cccc}1& v_j& ...& v_j^Q\end{array}\right]}^T,$ j＝1，2，...，J是一个大小为(Q+1)×1的实矢量，Q是多项式的阶数，v的大小为(Q+1)×J；

(3b3)按下式计算不同角度的扩展导向矢量在子空间每一维矢量的投影系数a_p(v_j)，

p＝1，...，P，j＝1，...，J，

其中u_p是低秩子空间U_P的第p维列矢量， $s_1\left(v_j\right)=e^{-i2\mathrm{\π}(2N-2)v_j}{\left[\begin{array}{cccc}1& e^{i2\mathrm{\π}{v}_j}& ...& e^{i2\mathrm{\π}(4N-4)v_j}\end{array}\right]}^T,$ j＝1，2，...，J是大小为(4N-3)×1的扩展导向矢量；

把不同维的投影系数按不同角度排成一列，得到大小为J×1投影系数矢量a_p，a_p＝[a_p(v₁) a_p(v₂) … a_p(v_J)]^T，p＝1，...，P；

(3b4)按下式计算子空间每一维矢量u_p的投影系数矢量a_p对应的多项式拟合系数矢量b_p：

b_p＝(VV^T)^-1Va_p，p＝1，...，P，

其中b_p的大小是(Q+1)×1，上标-1表示矩阵求逆，

把所有维的多项式拟合系数矢量排成矩阵，并求转置，得到大小为P×(Q+1)系数矩阵B：

B＝[b₁ b₂ … b_P]^T。

(3c)用步骤(3a2)求出的低秩子空间U_P、步骤(3b4)求出的系数矩阵B和对角阵D相乘，得到把扩展导向矢量转化为多项式的变换矩阵：

F＝U_PDB

其中D是大小为P×P的对角阵，P是低秩子空间U_P的维数，D＝diag([ε₁ ε₂ … ε_P])，diag表示对角阵，D的对角线元素ε_p满足

i为虚数单位。

步骤4：把实际波束宽度等分成4部分，每一部分的波束宽度小于虚拟波束宽度，运用步骤2求出的含有目标角度信息的复数据列矢量c、步骤3求出的变换矩阵F和与实际波束宽度等分有关的对角阵Λ_m，把雷达系统归一化空间频率v的似然函数J(v)的导数J(v)′转化为4组低阶多项式J_m(v)′＝β(v)^HF^HΛ_mc，其中m是波束宽度等分序号，m的取值是m＝1，2，3，4，

其中上标′表示求导，β(v)＝[1 v … v^Q]^T为多项式矢量，Q为多项式阶数，上标T表示转置，上标H表示共轭转置，

diag表示对角阵，s₁(v)＝e^-i2π(2N-2)v[1 e^i2πv … e^i2π(4N-4)v]^T是大小为(4N-3)×1的扩展导向矢量，i是虚数单位，v₀是视线方向归一化空间频率，

是实际波束宽度，N是阵元数，归一化空间频率v的取值范围是

步骤5：用牛顿法求上述4组低阶多项式的实根，目标的归一化空间频率可能存在这些实根中。

步骤6：把实际波束宽度的边缘点和步骤5中求出的实根代入归一化空间频率似然函数

其中s(v)＝[1 e^i2πv … e^i2π(N-1)v]^T是实际波束宽度内的导向矢量，计算出似然函数的最大值，最大值对应的归一化空间频率就为目标的归一化空间频率v_t的最大似然估计，用归一化空间频率和角度的转换关系，即可求出目标的角度：

$\mathrm{\θ}=\arcsin \left(\frac{v_t\mathrm{\λ}}{d}\right)$

其中θ为目标偏离阵列法线方向的角度，λ为雷达波长，d为雷达天线阵元间距。本发明是在主瓣和旁瓣均存在干扰背景下实现的，尤其对于大阵列雷达系统提供了一种快速测角新方法。

本发明通过对回波数据进行构造含有目标角度信息的复数据列矢量和把虚拟波束宽度的任意一个扩展导向矢量分解为变换矩阵和多项式矢量乘积的处理，使对目标归一化空间频率的似然函数直接求导成为可能，并且把导数转化为多项式，又通过把实际波束宽度等分成4部分，把所述导数用4组与阵元数无关的低阶多项式表示，降低了计算复杂度，保证了数值稳定性，而且适用于大阵列。

实施例2

主瓣和旁瓣干扰背景下雷达目标角度快速估计方法同实施例1，其中步骤2中构造含有目标角度信息的复数据列矢量c的方法更具体的步骤是：

(2a)利用所有距离单元的回波数据估计干扰和噪声的协方差矩阵R，并对协方差矩阵求逆，得到干扰和噪声的协方差矩阵的逆矩阵：

$R^{-1}={\left(\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{x}_l{x}_l^H\right)}^{-1}$

其中R^-1为干扰和噪声的协方差矩阵的逆矩阵，L为接收的雷达回波数据的距离单元数，L是由具体的雷达系统确定的，x_l是回波数据列矢量，大小为N×1，上标H表示共轭转置，上标-1表示矩阵求逆；

(2b)利用干扰和噪声的协方差矩阵的逆矩阵R^-1和实际波束宽度的主波束中心对应的导向矢量 $s\left(v_c\right)={\left[\begin{array}{cccc}1& e^{i2\mathrm{\π}{v}_c}& ...& e^{i2\mathrm{\π}(N-1)v_c}\end{array}\right]}^T$ 相乘，其中i为虚数单位，v_c为主波束中心的归一化空间频率，上标T表示转置，计算出自适应权值w＝R^-1s(v_c)，用自适应权值对雷达所有距离单元回波数据进行滤波，并用单元平均恒虚警方法确定出目标所在的距离单元，选取目标所在距离单元的雷达回波数据作为用来估计雷达目标角度的数据，为了方便说明，对目标所在距离单元的雷达回波数据建模，即：

x_t＝a_ts(v_t)+n

其中x_t为目标所在距离单元的雷达回波数据列矢量，a_t为目标复幅度， $s\left(v_t\right)={\left[\begin{array}{cccc}1& e^{i2\mathrm{\π}{v}_t}& ...& e^{i2\mathrm{\π}(N-1)v_t}\end{array}\right]}^T$ 是与目标归一化空间频率v_t有关的目标导向矢量，上标T表示转置，n是干扰加噪声矢量；

(2c)利用干扰和噪声的协方差矩阵的逆矩阵R^-1和目标所在距离单元的雷达回波数x_t据构造白化数据矩阵A，利用白化数据矩阵的对角元素构造列矢量1；

白化数据矩阵A按下式构造：

$A=R^{-1}{x}_t{x}_t^H{R}^{-1}$

利用白化数据矩阵A的对角线元素，按下两式构造列矢量1：

$a_k=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}A(N-k+n,n),$ k＝1，2，...，N

$a_1={\left[\begin{array}{ccccccc}{a}_1& \·\·\·& a_{N-1}& a_N& a_{N-1}^{*}& \·\·\·& a_1^{*}\end{array}\right]}^T$

其中a₁为大小为(2N-1)×1的列矢量1，上标*表示共轭，上标T表示转置，a_k是构造a₁的基本元素，k的取值范围为1，2，…，N，N为阵元数；

(2d)利用干扰和噪声的协方差矩阵的逆矩阵R^-1，按下两式构造列矢量2：

$r_k=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{R}^{-1}(N-k+n,n),$ k＝1，2，...，N

$a_2={\left[\begin{array}{ccccccc}{r}_1& \·\·\·& r_{N-1}& r_N& r_{N-1}^{*}& \·\·\·& r_1^{*}\end{array}\right]}^T$

其中a₂为大小为(2N-1)×1的列矢量2，r_k是构造a₂的基本元素，k的取值范围为1，2，…，N；

(2e)用列矢量1、列矢量2和对角阵D₁按下式求得一个新矩阵：

$C=D_1^H{a}_1{a}_2^H+a_1{a}_2^H{D}_1$

其中C是大小为(2N-1)×(2N-1)的新矩阵，D₁是对角线元素是首项为-i2π(N-1)、公差是i2π的等差数列的对角阵，即D₁＝diag([-i2π(N-1) … 0 … i2π(N-1)])，diag表示对角阵，i为虚数单位。

(2f)利用新矩阵C，按下两式构造出含有目标角度信息的复数据列矢量c：

$c_k=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}C(N-k+n,n),$ k＝1，2，...，2N-1

$c={\left[\begin{array}{ccccccc}{c}_1& \·\·\·& c_{2N-2}& c_{2N-1}& c_{2N-2}^{*}& \·\·\·& c_1^{*}\end{array}\right]}^T$

其中c为大小为(4N-3)×1的含有目标角度信息的复数据列矢量，c_k是构造c的基本元素，k的取值范围为1，2，…，2N-1。

本发明通过对回波数据进行构造含有目标角度信息的复数据列矢量，初步获得了含有目标角度信息的数据，为下一步通过变换矩阵把似然函数的导数变换为多项式做准备。

实施例3

主瓣和旁瓣干扰背景下雷达目标角度快速估计方法同实施例1-2

本例通过仿真实验给出本发明中的子空间拟合方法对于不同的阵元数，子空间拟合误差与子空间维数的关系，主要涉及步骤(3a)的过程。

1.实验条件

仿真在MATLAB7.0软件下进行，仿真参数如下：视线方向v₀为法线方向，阵元数分别为10、20和30，对应的波束宽度分别为

和

2.实验结果

图2给出了对于不同的阵元数，子空间拟合误差随着子空间维数变化的曲线图。

从图2中可以看出，1)随着子空间维数的增加，子空间拟合误差快速减小，当子空间维数为6时，子空间拟合误差已经非常小，表明导向矢量矩阵可以由6个基矢量很好地表示；2)子空间拟合误差对于不同的阵元数是基本相同的，表明张成主波束空间的基矢量数量和阵元数是独立的，这对于减少计算复杂度是非常有用的性质。

实施例4

主瓣和旁瓣干扰背景下雷达目标角度快速估计方法同实施例1-2

本例通过仿真实验给出本发明中的多项式拟合方法对于不同的阵元数，多项式拟合误差与多项式阶数的关系，主要涉及步骤(3b)的过程。

1.实验条件

仿真在MATLAB7.0软件下进行，仿真参数如下：视线方向v₀为法线方向，阵元数分别为10、20和30，对应的波束宽度分别为

和

2.实验结果

图3给出了对于不同的阵元数，多项式拟合误差随着多项式阶数变化的曲线图。

从图3中可以看出：1)随着多项式阶数的增加，多项式拟合误差快速减小，当多项式阶数小于10时，多项式拟合误差已经非常小，表明扩展导向矢量在低秩子空间的投影系数矢量可以用低阶多项式很好地近似；2)多项式拟合误差对于不同的阵元数是基本相同的，表明多项式的阶数和阵元数是独立的。对于实际测角中采用的大阵列，传统的求根最大似然方法计算复杂度与阵列的阵元数有关，当阵列阵元数较小时，测角计算复杂度和本发明的方法相差不多，但是随着雷达阵列阵元数增加，传统求根方法的计算复杂度随之增大，计算过程长，效率低，不利于工程化实际运用。本发明提出的求根最大似然方法中，多项式的阶数和阵元数无关，测角计算复杂度也和阵元数无关，因此，本发明有效地降低了计算复杂度，无论阵列的阵元数有多大，本发明的方法都能快速地测出目标角度。

实施例5

主瓣和旁瓣干扰背景下雷达目标角度快速估计方法同实施例1-2

本例通过单次仿真实验比较本发明提出的求根最大似然方法和前人提出的网格搜索最大似然方法、传统的求根最大似然方法的性能。

1.实验条件

雷达天线为一个均匀线阵，阵元数为20，阵元间距为半波长，视线方向指向阵面法向，干噪比为40dB的主瓣窄带干扰从偏离法线方向-1°入射，干噪比定义为

其中N是阵元数，

和

分别是阵元级干扰和噪声功率。另外，4个干噪比为30dB的旁瓣窄带干扰分别从偏离法线方向[-27°-9° 10°30°]入射。一个信噪比为25dB的目标从偏离法线方向2°入射，对于给定目标的信噪比定义为a_t为目标复幅度。采用只包含干扰加噪声独立同分布的200个距离单元的回波数据估计干扰和噪声的协方差矩阵。子空间的维数P为6，多项式的阶数Q为7，用于网格搜索最大似然估计的角度样本数为500。

2.实验结果

图4是本发明仿真的在实际波束宽度内网格搜索最大似然方法似然函数值的曲线图及传统求根最大似然方法似然函数的极值点和本发明的求根最大似然方法似然函数的极值点图。

在图4中，传统求根最大似然方法的极值点用o标出，本发明提出的求根最大似然方法的极值点用*标出，两种方法求出的极值点是相同的，并且这些极值点都落在最优的网格搜索最大似然方法似然函数值的曲线上，表明归一化空间频率的似然函数的极值点可以由本发明的求根最大似然估计有效地确定，似然函数取得最大值的角度即为目标角度，从而证明了本发明和网格搜素最大似然方法、传统求根最大似然的方法性能相同，但是本发明所需的计算复杂度较低，在相同的测角性能下，本发明更高效。

实施例6

主瓣和旁瓣干扰背景下雷达目标角度快速估计方法同实施例1-2

本例通过测量21个目标角度，每个目标角度测量值统计300次的仿真实验，对本发明提出的求根最大似然方法和前人提出的网格搜索最大似然方法、传统的求根最大似然方法的性能进行比较。

1.实验条件

雷达天线为一个均匀线阵，阵元数为20，阵元间距为半波长，视线方向指向阵面法向，干噪比为40dB的主瓣窄带干扰从偏离法线方向-1°入射，干噪比定义为

其中N是阵元数，和分别是阵元级干扰和噪声功率。另外，4个干噪比为30dB的旁瓣窄带干扰分别从偏离法线方向[-27°-9° 10°30°]入射。信噪比为25dB的目标分别从偏离法线方向的以下21个角度依次[-2.866° -2.5792° -2.2924° -2.0058° -1.7191° -1.4325° -1.146° -0.8595° -0.5730-0.2865° 0° 0.2865° 0.573° 0.8595° 1.146° 1.4325° 1.7191° 2.0058° 2.2924°2.5792° 2.866°]入射，对于给定目标的信噪比定义为a_t为目标复幅度。采用只包含干扰加噪声独立同分布的200个距离单元的回波数据估计干扰和噪声的协方差矩阵。子空间的维数P为6，多项式的阶数Q为7，用于网格搜索最大似然估计的角度样本数为500。

2.实验结果

图5是本发明与网格搜索最大似然方法、传统求根最大似然方法针对目标角度和均方根误差的曲线比较图，横坐标是目标角度，纵坐标是均方根误差。

从图5中看出：通过300次的统计平均，本发明提出的求根最大似然方法角度的均方根误差和网格搜索最大似然方法、传统的求根最大似然方法角度的均方根误差在所有估计出的目标角度中几乎没有差别。因此，本发明提出的求根最大似然方法角度测量性能与网格搜索最大似然方法、传统的求根最大似然方法性能基本一致，并且通过300次的统计试验，试验结果相同，说明本发明提出的求根最大似然方法测角稳健，有好的数值稳定性。

综上所述，本发明的主瓣和旁瓣干扰背景下雷达目标角度快速估计方法，其步骤简述为：(1)接收雷达回波数据；(2)由雷达回波数据得到含有目标角度信息的复数据列矢量；(3)计算把扩展导向矢量变换为多项式的变换矩阵；(4)把雷达系统归一化空间频率的似然函数的导数转化为4组低阶多项式；(5)用牛顿法求4组低阶多项式的实根；(6)在实际波束宽度的边缘点和4组多项式的根中，求出目标归一化空间频率，得到目标角度。本发明通过对回波数据进行构造含有目标角度信息的复数据列矢量和把虚拟波束内任意一个扩展导向矢量分解为变换矩阵和多项式矢量乘积形式的处理，使对目标归一化空间频率的似然函数直接求导成为可能，并且通过把实际波束宽度等分成4部分，把目标归一化空间频率的似然函数求导用4组低阶实多项式表示，从而似然函数的极点由实际波束宽度边缘点和低阶实多项式的根确定，克服了现有技术计算复杂度高、数值稳定性差和不适用于大阵列的不足，具有计算复杂度低，数值稳定性好，特别适用于大阵列的优点。

Claims (4)

1.一种主瓣和旁瓣干扰背景下雷达目标角度快速估计方法，其特征在于包括如下步骤：

步骤1：雷达通过天线接收回波数据，该数据按阵元通道和距离单元二维矩阵的格式存储，回波数据中包括主瓣干扰和旁瓣干扰、目标信号和接收机噪声；

步骤2：用雷达回波数据估计干扰和噪声的协方差矩阵，用自适应滤波、单元平均恒虚警检测方法确定目标所在的距离单元，用干扰和噪声的协方差矩阵的逆矩阵和目标所在距离单元的雷达回波数据构造白化数据矩阵A，用该白化数据矩阵信息和所述协方差矩阵的逆矩阵信息构造两个列矢量，利用两个列矢量信息得到含有目标角度信息的复数据列矢量c；

步骤3：计算把扩展导向矢量变换为多项式的变换矩阵F；

(3a)求雷达虚拟波束宽度内扩展导向矢量的低秩子空间U_P；

(3b)在虚拟波束宽度内，把不同角度的扩展导向矢量在子空间的投影系数矢量用多项式进行拟合，得到系数矩阵B；

(3c)用低秩子空间U_P、系数矩阵B和对角阵D相乘，得到把扩展导向矢量转化为多项式的变换矩阵：

F＝U_PDB

其中D是大小为P×P的对角阵，P是低秩子空间U_P的维数，D＝diag([ε₁ ε₂ … ε_P])，diag表示对角阵，D的对角线元素ε_p满足

i为虚数单位；

步骤4：把实际波束宽度等分成4部分，每一部分的波束宽度小于虚拟波束宽度，运用含有目标角度信息的复数据列矢量c、变换矩阵F和与实际波束宽度等分有关的对角阵Λ_m，把雷达系统归一化空间频率v的似然函数J(v)的导数J(v)′转化为4组低阶多项式J_m(v)′＝β(v)^HF^HΛ_mc，m＝1，2，3，4，

其中上标′表示求导，β(v)＝[1 v … v^Q]^T为多项式矢量，Q为多项式阶数，上标T表示转置，上标H表示共轭转置，

diag表示对角阵，s₁(v)＝e^-i2π(2N-2)v[1 e^i2πv … e^i2π(4N-4)v]^T是大小为4N-3的扩展导向矢量，i是虚数单位，v₀是视线方向归一化空间频率，

是实际波束宽度，N是阵元数，归一化空间频率v的取值范围是 $v\∈[v_0-\frac{1}{2}{\mathrm{\Δ}}_N,v_0+\frac{1}{2}{\mathrm{\Δ}}_N];$

步骤5：用牛顿法求4组低阶多项式的实根；

步骤6：把实际波束宽度的边缘点和4组低阶多项式的实根代入归一化空间频率似然函数J(v)，计算出似然函数的最大值，最大值对应的归一化空间频率就为目标的归一化空间频率的最大似然估计，用归一化空间频率和角度的转换关系，求出目标的角度。

2.根据权利要求1所述的主瓣和旁瓣干扰背景下雷达目标角度快速估计方法，其特征在于，步骤2中构造含有目标角度信息的复数据列矢量c的方法包括：

(2a)利用干扰和噪声的协方差矩阵的逆矩阵R^-1和目标所在距离单元的雷达回波数x_t据构造白化数据矩阵A，利用白化数据矩阵的对角元素构造列矢量1；

白化数据矩阵A按下式构造：

$A=R^{-1}{x}_t{x}_t^H{R}^{-1}$

其中上标-1表示矩阵求逆，上标H表示共轭转置，

利用白化数据矩阵A的对角线元素，按下两式构造列矢量1：

$a_k=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}A(N-k+n,n),$ k＝1，2，...，N

$a_1={\left[\begin{array}{ccccccc}{a}_1& \·\·\·& a_{N-1}& a_N& a_{N-1}^{*}& \·\·\·& a_1^{*}\end{array}\right]}^T$

其中a₁为大小为(2N-1)×1的列矢量1，上标*表示共轭，上标T表示转置，a_k是构造a₁的基本元素，k的取值范围为1，2，…，N，N为阵元数；

(2b)利用干扰和噪声的协方差矩阵的逆矩阵R^-1，按下两式构造列矢量2：

$r_k=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}{R}^{-1}(N-k+n,n),$ k＝1，2，...，N

$a_2={\left[\begin{array}{ccccccc}{r}_1& \·\·\·& r_{N-1}& r_N& r_{N-1}^{*}& \·\·\·& r_1^{*}\end{array}\right]}^T$

其中a₂为大小为(2N-1)×1的列矢量2，r_k是构造a₂的基本元素，k的取值范围为1，2，…，N；

(2c)用列矢量1、列矢量2和对角阵D₁按下式求得一个新矩阵：

$C=D_1^H{a}_1{a}_2^H+a_1{a}_2^H{D}_1$

其中C是大小为(2N-1)×(2N-1)的新矩阵，D₁是对角线元素是首项为-i2π(N-1)、公差是i2π的等差数列的对角阵，即D₁＝diag([-i2π(N-1) … 0 … i2π(N-1)])，diag表示对角阵，i为虚数单位。

(2d)利用新矩阵C，按下两式构造出含有目标角度信息的复数据列矢量c：

$c_k=\underset{n=1}{\overset{k}{\mathrm{\Σ}}}C(N-k+n,n),$ k＝1，2，...，2N-1

$c={\left[\begin{array}{ccccccc}{c}_1& \·\·\·& c_{2N-2}& c_{2N-1}& c_{2N-2}^{*}& \·\·\·& c_1^{*}\end{array}\right]}^T$

其中c为大小为(4N-3)×1的含有目标角度信息的复数据列矢量，c_k是构造c的基本元素，k的取值范围为1，2，…，2N-1。

3.根据权利要求2所述的主瓣和旁瓣干扰背景下雷达目标角度快速估计方法，其特征在于，步骤(3a)中得到雷达虚拟波束宽度内扩展导向矢量的低秩子空间U_P的过程包括：

(3a1)雷达天线阵元数为N，虚拟波束宽度为

在虚拟波束宽度内构造sinc矩阵，记为G，它的每一个元素按下式计算：

$g_{\mathrm{xy}}\left(\mathrm{\Δ}\right)=\mathrm{\Δ}{e}^{i2\mathrm{\π}(x-y)v_0}\sin c\left(\mathrm{\Δ}(x-y)\right),$ x＝1，2，…，4N-3；y＝1，2，…，4N-3；

(3a2)对sinc矩阵G做特征分解，选取大的特征值对应的特征向量构成U_P，选取的特征向量的个数为6，即子空间的维数P为6。

4.根据权利要求3所述的主瓣和旁瓣干扰背景下雷达目标角度快速估计方法，其特征在于，步骤(3b)中把虚拟波束宽度内，不同角度的扩展导向矢量在低秩子空间的投影系数矢量用多项式进行拟合，得到系数矩阵B的过程包括：

(3b1)在虚拟波束宽度内取不同角度的样本点v₁，v₂，...，v_J，其中J是选取角度样本点的个数；

(3b2)按下式构造不同角度对应的多项式矩阵v：

V＝[β(v₁) β(v₂) ... β(v_J)]

其中 $\mathrm{\β}\left(v_j\right)={\left[\begin{array}{cccc}1& v_j& ...& v_j^Q\end{array}\right]}^T,$ j＝1，2，...，J是大小为(Q+1)×1的实矢量，Q是多项式的阶数，v的大小为(Q+1)×J；

(3b3)按下式计算不同角度的扩展导向矢量在子空间每一维矢量的投影系数a_p(v_j)，

p＝1，...，P，j＝1，...，J，

其中u_p是低秩子空间U_P的第p维列矢量， $s_1\left(v_j\right)=e^{-i2\mathrm{\π}(2N-2)v_j}{\left[\begin{array}{cccc}1& e^{i2\mathrm{\π}{v}_j}& ...& e^{i2\mathrm{\π}(4N-4)v_j}\end{array}\right]}^T,$ j＝1，2，...，J是大小为(4N-3)×1的扩展导向矢量；

把不同维的投影系数按不同角度排成一列，得到大小为J×1投影系数矢量a_p，a_p＝[a_p(v₁) a_p(v₂) … a_p(v_J)]^T，p＝1，...，P；

(3b4)按下式计算子空间每一维矢量u_p的投影系数矢量a_p对应的多项式拟合系数矢量b_p：

b_p＝(VV^T)^-1Va_p，p＝1，...，P，

b_p的大小是(Q+1)×1，上标-1表示矩阵求逆，

把所有维的多项式拟合系数矢量排成矩阵，并求转置，得到大小为P×(Q+1)系数矩阵B：

B＝[b₁ b₂ … b_P]^T。

Images