CN104181513A - 一种雷达天线阵元位置的校正方法 - Google Patents

Abstract

本发明公开了一种雷达天线阵元位置的校正方法，涉及阵列信号处理领域，其步骤为：步骤1，构造回波信号模型矩阵并求其自相关矩阵,然后设定N个阵元的位置坐标；步骤2，估计出目标的到达角，并求出到达角相关矩阵的特征值对应的特征向量构成的矩阵；步骤3，估计阵元的扰动矩阵；步骤4，计算出阵元的位置扰动矩阵；步骤5，计算出阵元的位置坐标。本发明主要解决了阵元位置扰动难以估计的问题。本发明可以比较精确地估计出阵元的位置扰动，进而估计出阵元的位置。

Description

一种雷达天线阵元位置的校正方法

技术领域

本发明属于阵列信号处理领域，涉及一种雷达天线阵元位置的校正方法。

背景技术

为了估计出信号的到达角DOA(Direction Of Arrival)，然后估计出阵元的位置偏差,进而估计出阵元的位置，国内外进行了很多的研究。现在一种应用广泛的方法是多重信号分类MUSIC(Multiple Signal Classification)算法,Schmidt R O等人在1979年提出了MUSIC算法。MUSIC算法的基本思想是将任意阵列的输出数据的协方差矩阵进行特征值分解，从而得到与信号分量相对应的信号子空间和与信号分量正交的噪声子空间，然后利用这两个子空间的正交性来估计信号的参数(入射方向、极化信息及信号强度等)，从而可以估计出信号的DOA。

经过仿真发现，MUSIC算法虽然明显提高了信号的分辨率，但是仍然会有偏差，当要求估计出精确地DOA时MUSIC算法仍然有它自己的局限性。同时，在实际的工程应用中，由于各种误差不可避免，实际的阵列流型往往会出现一定的偏差或扰动，此时，通常的高分辨空间谱(如：MUSIC算法估计出的高分辨谱)估计算法的性能会严重恶化，甚至失效，因此，阵列误差一直是高分辨空间谱估计技术受到限制的一个重要原因。于是，人们开始研究阵列误差的校正。早期的阵列校正正是通过对阵列流型直接进行离散测量、内插、存储来实现的，但这些方法实现代价较大而且效果并不明显。

因此，20世纪90年代以后人们通过对阵列扰动进行建模，将阵列误差校正逐渐转化为一个参数估计问题。参数类的阵列校正方法通常可以分为有源校正类和自校正类。有源校正通过在空间设置方位精确已知的辅助信源对阵列扰动参数进行联合估计，而自校正类方法通常根据某种优化函数对空间信源的方位与阵列的扰动参数进行联合估计。这两种算法各有优缺点：对于有源校正而言，不需要对信号源方位进行估计，所以其运算量较小，因此实际中被采纳的比较多，但这类校正算法对辅助信号源有较高的精确方位信息要求；而自校正算法可以不需要方位已知的辅助信源，而且可以在线完成实际方位估计，所以校正的精度比较高，但由于误差参数与方位参数之间的耦合和某些病态的阵列结构，参数估计的唯一辨识通常无法保证。因此在使用哪种方法来对阵列位置误差进行校正时要考虑实际情况，本专利采用了自校正算法。在自校正算法中需要用到迭代，设代价函数为J，在每一方位和频率迭代估计的每一步中都会减小，同时由于J≥0，所以该最优化过程可以保持收敛到一个局部最优点，但不一定是全局最优点。于是，自校正算法也存在自己的局限性：

(1)对于均匀线阵，由于其导向矢量的范德蒙特性，方位估计与相位误差估计存在模糊性；

(2)阵元数小于4时算法失效，而且阵元数大于4时，对于某些特殊的阵列结构和方位组合，算法的解也可能不唯一。

因此无论是采用多重信号分类(MUSIC)算法还是自校正算法，都无法精确估计出阵元的位置，都会产生一定的误差。由于这些缺点，本专利提出了下面的方法来估计阵元的位置，对实际很有指导意义。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种雷达天线阵元位置的校正方法，能够估计出阵元的位置偏差，并进一步估计出阵元的位置。

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

一种雷达天线阵元位置的校正方法，其特征在于，包括以下步骤：

步骤1，雷达天线接收原始的回波信号，建立该回波信号的模型矩阵Z；利用回波信号的模型矩阵Z求出回波信号的自相关矩阵R；并且设定雷达天线的N个阵元的位置矩阵(X',Y')，其中(X',Y')＝[(X'₁,Y'₁)(X'₂,Y'₂)...(X'_N,Y'_N)]^T表示N个阵元的位置矩阵；

步骤2，利用回波信号的模型矩阵Z和回波信号的自相关矩阵R估计出目标的到达角θ；构造到达角相关矩阵Q(θ)，并求出到达角相关矩阵Q(θ)的特征值对应的特征向量构成矩阵V；

步骤3，根据目标的到达角θ和到达角相关矩阵Q(θ)的特征值对应的特征向量构成矩阵V，估计出N个阵元的扰动矩阵Γ₁(θ)，Γ₁(θ)＝[(Γ₁(θ))₁,(Γ₁(θ))₂,...,(Γ₁(θ))_N)]^T；

步骤4，利用估计出的N个阵元的扰动矩阵Γ₁(θ)，估计出N个阵元的位置偏差矩阵(ΔX₁,ΔY₁)，(ΔX₁,ΔY₁)＝[(ΔX₁,ΔY₁)₁,(ΔX₁,ΔY₁)₁,...,((ΔX₁,ΔY₁)_N)]^T；

步骤5，利用估计出的N个阵元的位置偏差矩阵(ΔX₁,ΔY₁)，得到估计出的阵元位置(X₁,Y₁)，估计出的阵元位置(X₁,Y₁)＝(X',Y')+(ΔX₁,ΔY₁)。

上述技术方案的特点和进一步改进在于：

(1)步骤1包括以下子步骤：

1a)构建目标的导向矢量矩阵A(θ)：

$A\left(\mathrm{\θ}\right)=\exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}\×(0:N-1)\×d\×\sin (\mathrm{\θ}\×\frac{\mathrm{\π}}{180}))---\left(1\right)$

其中,N表示雷达天线的阵元数，λ表示雷达工作的波长，d表示阵元的间隔距离,θ表示目标的到达角，θ＝[θ₁,θ₂,...θ_l,...,θ_K],θ_l表示第l个目标的到达角，K表示目标的个数；

1b)利用目标的导向矢量矩阵A(θ)，构建雷达天线的N个阵元的初始扰动矩阵Γ(θ)：

其中，Γ(θ)为N×1维矩阵，N为雷达天线的阵元数，A(θ)表示目标的导向矢量矩阵，表示阵元相位的随机扰动矩阵，θ表示目标的到达角；

1c)利用目标的导向矢量矩阵A(θ)和初始扰动矩阵Γ(θ)构造回波信号的模型矩阵Z；

Z＝Γ(θ)A(θ)+E      (3)

其中，A(θ)表示目标的导向矢量矩阵，Γ(θ)表示N个阵元的初始扰动矩阵，E表示噪声矩阵，θ表示目标的到达角；

利用回波信号的模型矩阵Z求出回波信号的自相关矩阵R：

$R=\frac{1}{N}{\mathrm{ZZ}}^H---\left(4\right)$

其中，Z表示回波信号的模型矩阵，N表示雷达天线的阵元数，[·]^H表示共轭转置；

1d)设定雷达天线的N个阵元的位置矩阵(X',Y')，其中，(X',Y')＝[(X₁',Y₁'),(X₂',Y₂'),...,(X_N',Y_N')]^T。

(2)步骤2包括以下子步骤：

2a)对目标的自相关矩阵R进行特征值分解，得到目标子空间U_s和噪声子空间U_N；

2b)利用目标的导向矢量矩阵A(θ)和阵元的初始扰动矩阵Γ(θ)构造目标辅助矩阵w(θ)；

w(θ)＝a_m(θ)Γ(θ)      (5)

其中，a_m(θ)由A(θ)的第m列构成，m＝1,2,3...,K，K表示目标的个数，θ表示目标的到达角；

2c)根据目标辅助矩阵w(θ)与噪声子空间U_N的正交性，得到下式：

w^H(θ)U_NU^H _Nw(θ)＝0      (6)

其中，U_N表示噪声子空间，N表示雷达天线阵元数，[·]^H表示共轭转置，θ表示目标的到达角；

将目标辅助矩阵w(θ)公式(5)代入公式(6)得：

Γ^H(θ)a_m ^H(θ)U_NU^H _Na_m(θ)Γ(θ)＝0      (7)

由公式(7)得到达角相关矩阵Q(θ)：

Q(θ)＝a_m ^H(θ)U_NU^H _Na_m(θ)      (8)

2d)将到达角相关矩阵Q(θ)代入以下公式(9)求出目标的到达角θ：

$\mathrm{\θ}=\arg \underset{\mathrm{\θ}}{\max }\frac{1}{{\mathrm{\λ}}_{\min }\left[Q\left(\mathrm{\θ}\right)\right]}$ 或 $\mathrm{\θ}=\arg \underset{\mathrm{\θ}}{\max }\frac{1}{\det \left[Q\left(\mathrm{\θ}\right)\right]}---\left(9\right)$

其中，arg[·]为求解最优化，max表示求最大值，λ_min[·]为求矩阵的最小特征值，[·]^H表示矩阵的共轭转置，det[·]为求矩阵的行列式，θ表示目标的到达角；

2e)对到达角相关矩阵Q(θ)进行奇异值分解，得到达角相关矩阵Q(θ)的特征值对应的特征向量构成矩阵V。

(3)步骤3包括以下子步骤：

3a)将求解阵元的扰动矩阵Γ₁(θ)转换为如下的优化模型公式(10)：

${\mathrm{\Γ}}_1\left(\mathrm{\θ}\right)=\underset{S.D}{\min }{\left|\right|Z-[A\left(\mathrm{\θ}\right)+D]S\left|\right|}_F^2+\mathrm{\μ}{\left|\right|S\left|\right|}_{\∞,0}---\left(10\right)$

其中，||·||_F为Frobenius范数运算符，Z表示回波信号的模型矩阵，A(θ)是目标的导向矢量，S为目标的稀疏矩阵，D为阵列流行矩阵；||·||_∞，0表示混合范数，μ>0为正则化参数，θ表示目标的到达角；

3b)设定阵列流型矩阵D等于Q(θ)的特征值对应的特征向量构成矩阵V，即：

D＝V      (11)

3c)将优化模型公式(10)转换为以下模型公式(12)：

$\underset{S}{\min }{\left|\right|Z-[A\left(\mathrm{\θ}\right)+D]S\left|\right|}_F^2---\left(12\right)$

其中，Z表示回波信号的模型矩阵，A(θ)表示目标的导向矢量矩阵，S表示Q(θ)的非零行构成的稀疏矩阵，D表示阵列流行矩阵，D等于到达角相关矩阵Q(θ)的特征值对应的特征向量构成矩阵V，θ表示目标的到达角；

3d)设定Φ＝A(θ)D，然后设定Φ的支撑矩阵Ω，支撑矩阵Ω中包含了Φ中不为零的列；

3e)β是支撑矩阵Ω中的一个元素，通过求解式(13)得到β的第i次迭代元素β_i：

${\mathrm{\β}}_i=\arg \min \left[\mathrm{\Σ}\right|({\mathrm{\Φ}}_{{\mathrm{\Ω}}_{i-1}}^H{P}_{i-1}-{\left|\right|{\mathrm{\Φ}}_{{\mathrm{\Ω}}_{i-1}}^H{P}_{i-1}\left|\right|}_{\∞})\left|\right]---\left(13\right)$

其中，表示在第i-1次迭代中由支撑矩阵Ω_i-1中的元素对应于矩阵Φ中的列向量构成的矩阵，设定支撑矩阵Ω的初值Ω₀＝Φ；

构造P_i-1表示矩阵P第i-1次迭代矩阵，设定矩阵P的初值P₀＝Z，[·]^H表示共轭转置，arg表示求解最优化问题，i＝1,2,...,K，K表示目标数；

3f)利用第i次迭代元素β_i求解第i次迭代时的支撑矩阵Ω_i＝Ω_i-1∪β_i；

3g)令迭代次数i增加1，重复迭代3e)-3f)步，直到i等于K，得到Ω_K迭代结束；K表示目标的个数，并设定支撑矩阵Ω＝Ω_K，Ω_K＝Ω_K-1∪β_k共有K个元素；

3h)利用求出的支撑矩阵Ω，求出稀疏矩阵S：

定义矩阵S_Ω，S_Ω计算公式为其中，Φ_Ω表示由支撑矩阵Ω中的元素对应于矩阵Φ中的列向量构成的矩阵，[·]^H表示共轭转置操作，Z表示回波信号的模型矩阵，[·]^-1表示求矩阵的逆，A(θ)是目标的导向矢量，θ表示目标的到达角；

矩阵S_Ω组成目标的稀疏矩阵S的非零列，令S的其余列为零，得到稀疏矩阵S；

3k)将阵列流型矩阵D和稀疏矩阵S代入优化模型为(10)中，估计出N个阵元的扰动矩阵Γ₁(θ)：

${\mathrm{\Γ}}_1\left(\mathrm{\θ}\right)={\left|\right|Z-[A\left(\mathrm{\θ}\right)+D]S\left|\right|}_F^2+\mathrm{\μ}{\left|\right|S\left|\right|}_{\∞,0}---\left(14\right)$

其中，其中，为Frobenius范数运算符，Z表示回波信号的模型矩阵，A(θ)是目标的导向矢量，S表示稀疏矩阵，D表示阵列流行矩阵，D等于到达角相关矩阵Q(θ)的特征值对应的特征向量构成矩阵V；||·||_∞，0表示混合范数，μ>0为正则化参数，θ表示目标的到达角。

(4)步骤4具体包括：

将步骤3估计的阵元的扰动矩阵Γ₁(θ)代入以下公式(15)中，求出N个阵元的位置偏差矩阵(ΔX₁,ΔY₁)：

$\left\{\begin{array}{c}{P}_1\left(\mathrm{\θ}\right)=-\frac{\mathrm{\λ}}{2\mathrm{\π}}\mathrm{angle}\left({\mathrm{\Γ}}_1\left(\mathrm{\θ}\right)\right)\\ ({\mathrm{\ΔX}}_1,{\mathrm{\ΔY}}_1)=\left[\begin{array}{cc}{P}_1\left({\mathrm{\θ}}_1\right)& P_1\left({\mathrm{\θ}}_2\right),...,P_1\left({\mathrm{\θ}}_l\right),...,P_1\left({\mathrm{\θ}}_K\right)\end{array}\right]{\left({\left[\begin{array}{cccccc}{\sin \mathrm{\θ}}_1& {\sin \mathrm{\θ}}_2& & {\sin \mathrm{\θ}}_l& & {\sin \mathrm{\θ}}_K\\ & & ...& & ...& \\ {\cos \mathrm{\θ}}_1& {\cos \mathrm{\θ}}_2& & {\cos \mathrm{\θ}}_l& & {\cos \mathrm{\θ}}_K\end{array}\right]}^T\right)}^{-1}\end{array}\right.---\left(15\right)$

其中，λ表示雷达的波长，[]^-1表示求矩阵的逆，Γ₁(θ)表示估计出的阵元的扰动矩阵，angle为求向量的角度，θ表示目标的到达角，K表示目标的个数。

与现有技术相比，本发明具有突出的实质性特点和显著的进步。本发明与现有方法相比，具有以下优点：

现有估计阵元位置的方法都会有一定的误差，而本发明提出的方法可以比较精确地估计出阵元扰动，然后可以估计出阵元的位置偏差，从而估计出阵元的位置。

附图说明

下面结合附图和具体实施方式对本发明做进一步说明。

图1是本发明的实现流程图；

图2是分别利用现有技术中的MUSIC算法、DBF算法以及自校正算法估计出的目标的到达角的，图(a)是现有技术的三种方法计算得到的目标的到达角，图(b)是图(a)的局部放大图，横坐标表示目标的到达角单位为度，纵坐标为幅度，单位为dB。

图3是将实际的阵元位置偏差和用本发明提出的算法估计出的阵元的位置偏差进行仿真比较，横坐标表示阵元数，纵坐标表示阵元的位置偏差，单位为米。

图4是用实际的阵元位置偏差和用本发明提出的算法估计出的阵元的位置偏差分别对阵元位置进行校正的对比图，仿真图中横坐标表示校正后的阵元的水平位置，纵坐标表示校正后的阵元的高度。

具体实施方式

参照图1，说明本发明的一种雷达天线阵元位置的校正方法，包括以下步骤：

步骤1，雷达天线接收原始的回波信号，建立该回波信号的模型矩阵Z；利用回波信号的模型矩阵Z求出回波信号的自相关矩阵R；并且设定雷达天线的N个阵元的位置矩阵(X',Y')，其中(X',Y')＝[(X'₁,Y'₁)(X'₂,Y'₂)...(X'_N,Y'_N)]^T表示N个阵元的位置矩阵。

1a)构建目标的导向矢量矩阵A(θ)：

$A\left(\mathrm{\θ}\right)=\exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}\×(0:N-1)\×d\×\sin (\mathrm{\θ}\×\frac{\mathrm{\π}}{180}))---\left(1\right)$

其中,N表示雷达天线的阵元数，λ表示雷达工作的波长，d表示阵元的间隔距离,θ表示目标的到达角，θ＝[θ₁,θ₂,...θ_l,...,θ_K],θ_l表示第l个目标的到达角，K表示目标的个数；

1b)利用目标的导向矢量矩阵A(θ)，构建雷达天线的N个阵元的初始扰动矩阵Γ(θ)：

其中，Γ(θ)为N×1维矩阵，N为雷达天线的阵元数，A(θ)表示目标的导向矢量矩阵，表示阵元相位的随机扰动矩阵，θ表示目标的到达角；

1c)利用目标的导向矢量矩阵A(θ)和初始扰动矩阵Γ(θ)构造回波信号的模型矩阵Z；

Z＝Γ(θ)A(θ)+E      (3)

其中，A(θ)表示目标的导向矢量矩阵，Γ(θ)表示N个阵元的初始扰动矩阵，E表示噪声矩阵，θ表示目标的到达角；

利用回波信号的模型矩阵Z求出回波信号的自相关矩阵R：

$R=\frac{1}{N}{\mathrm{ZZ}}^H---\left(4\right)$

其中，Z表示回波信号的模型矩阵，N表示雷达天线的阵元数，[·]^H表示共轭转置；

1d)设定雷达天线的N个阵元的位置矩阵(X',Y')，其中，(X',Y')＝[(X₁',Y₁'),(X₂',Y₂'),...,(X_N',Y_N')]^T。

步骤2，利用回波信号的模型矩阵Z和回波信号的自相关矩阵R估计出目标的到达角θ；构造到达角相关矩阵Q(θ)，并求出到达角相关矩阵Q(θ)的特征值对应的特征向量构成矩阵V。

2a)对目标的自相关矩阵R进行特征值分解，得到目标子空间U_s和噪声子空间U_N；

2b)利用目标的导向矢量矩阵A(θ)和阵元的初始扰动矩阵Γ(θ)构造目标辅助矩阵w(θ)；

w(θ)＝a_m(θ)Γ(θ)      (5)

其中，a_m(θ)由A(θ)的第m列构成，m＝1,2,3...,K，K表示目标的个数，θ表示目标的到达角；

2c)根据目标辅助矩阵w(θ)与噪声子空间U_N的正交性，得到下式：

w^H(θ)U_NU^H _Nw(θ)＝0      (6)

其中，U_N表示噪声子空间，N表示雷达天线阵元数，[·]^H表示共轭转置，θ表示目标的到达角；

将公式(5)代入公式(6)得：

Γ^H(θ)a_m ^H(θ)U_NU^H _Na_m(θ)Γ(θ)＝0      (7)

由(7)式得到达角相关矩阵Q(θ)：

Q(θ)＝a_m ^H(θ)U_NU^H _Na_m(θ)      (8)

2d)将到达角相关矩阵Q(θ)代入公式(9)求出目标的到达角θ：

$\mathrm{\θ}=\arg \underset{\mathrm{\θ}}{\max }\frac{1}{{\mathrm{\λ}}_{\min }\left[Q\left(\mathrm{\θ}\right)\right]}$ 或 $\mathrm{\θ}=\arg \underset{\mathrm{\θ}}{\max }\frac{1}{\det \left[Q\left(\mathrm{\θ}\right)\right]}---\left(9\right)$

其中，arg[·]为求解最优化，max表示求最大值，λ_min[·]为求矩阵的最小特征值，[·]^H表示矩阵的共轭转置，det[·]为求矩阵的行列式，θ表示目标的到达角。

2e)对到达角相关矩阵Q(θ)进行奇异值分解，得到达角相关矩阵Q(θ)的特征值对应的特征向量构成矩阵V。

到达角相关矩阵Q(θ)的特征值对应的特征向量构成矩阵V的求取方法，例如：

i)目标的个数为K个，则θ＝[θ₁,θ₂,...θ_l,...,θ_K],θ_l表示第l个目标的到达角；

ii)到达角相关矩阵Q(θ)＝[Q(θ₁),Q(θ₂),...,Q(θ_K)]；

iii)对Q(θ)进行奇异值分解，得到特征值η＝[η₁,η₂,...,η_K]，特征值η＝[η₁,η₂,...,η_K]对

应的特征向量为V₁,V₂,...,V_K；

iv)用Q(θ)的特征值对应的特征向量V₁,V₂,...,V_K构成矩阵V＝[V₁,V₂,...,V_K]。

MATLAB中对到达角相关矩阵Q(θ)进行奇异值分解的方法为：

[V,η]＝eig(Q(θ)),其中eig(·)表示奇异值分解。

在步骤2中求取目标的到达角的方法即为自校正算法。

步骤3，根据目标的到达角θ和到达角相关矩阵Q(θ)的特征值对应的特征向量构成矩阵V，估计出N个阵元的扰动矩阵Γ₁(θ)，Γ₁(θ)＝[(Γ₁(θ))₁,(Γ₁(θ))₂,...,(Γ₁(θ))_N)]^T。

3a)将求解阵元的扰动矩阵Γ₁(θ)转换为如下的优化模型公式(10)：

${\mathrm{\Γ}}_1\left(\mathrm{\θ}\right)=\underset{S.D}{\min }{\left|\right|Z-[A\left(\mathrm{\θ}\right)+D]S\left|\right|}_F^2+\mathrm{\μ}{\left|\right|S\left|\right|}_{\∞,0}---\left(10\right)$

其中，||·||_F为Frobenius范数运算符，Z表示回波信号的模型矩阵，A(θ)是目标的导向矢量，S为目标的稀疏矩阵，D为阵列流行矩阵；||·||_∞，0表示混合范数，μ>0为正则化参数，θ表示目标的到达角；

3b)设定最优化模型(10)中的阵列流型矩阵D等于Q(θ)的特征值对应的特征向量构成矩阵V，即：

D＝V      (11)

以下子步骤3c)-3h)利用3b)中设定的阵列流型矩阵D，求解稀疏矩阵S；

3c)将优化模型公式(10)转换为以下模型公式(12)：

$\underset{S}{\min }{\left|\right|Z-[A\left(\mathrm{\θ}\right)+D]S\left|\right|}_F^2---\left(12\right)$

其中，Z表示回波信号的模型矩阵，A(θ)表示目标的导向矢量矩阵，S表示Q(θ)的非零行构成的稀疏矩阵，D等于到达角相关矩阵Q(θ)的特征值对应的特征向量构成矩阵V，θ表示目标的到达角；

3d)设定Φ＝A(θ)D，然后设定Φ的支撑矩阵Ω，支撑矩阵Ω中包含了Φ中不为零的列；

3e)β是支撑矩阵Ω中的一个元素，通过求解式(13)得到β的第i次迭代元素β_i：

${\mathrm{\β}}_i=\arg \min \left[\mathrm{\Σ}\right|({\mathrm{\Φ}}_{{\mathrm{\Ω}}_{i-1}}^H{P}_{i-1}-{\left|\right|{\mathrm{\Φ}}_{{\mathrm{\Ω}}_{i-1}}^H{P}_{i-1}\left|\right|}_{\∞})\left|\right]---\left(13\right)$

其中，表示在第i-1次迭代中由支撑矩阵Ω_i-1中的元素对应于矩阵Φ中的列向量构成的矩阵，设定支撑矩阵Ω的初值Ω₀＝Φ；

构造P_i-1表示矩阵P第i-1次迭代矩阵，设定矩阵P的初值P₀＝Z，[·]^H表示共轭转置，arg表示求解最优化问题，i＝1,2,...,K，K表示目标数；

3f)利用第i次迭代元素β_i求解第i次迭代时的支撑矩阵Ω_i＝Ω_i-1∪β_i；

3g)令迭代次数i增加1，重复迭代3e)-3f)步，直到i等于K，得到Ω_K迭代结束。K表示目标的个数，并设定支撑矩阵Ω＝Ω_K，Ω_K＝Ω_K-1∪β_K共有K个元素；

3h)利用求出的支撑矩阵Ω，求出稀疏矩阵S：

定义矩阵S_Ω，S_Ω计算公式为其中，Φ_Ω表示由支撑矩阵Ω中的元素对应于矩阵Φ中的列向量构成的矩阵，[·]^H表示共轭转置操作，Z表示回波信号的模型矩阵，[·]^-1表示求矩阵的逆，A(θ)是目标的导向矢量，θ表示目标的到达角；

矩阵S_Ω组成目标的稀疏矩阵S的非零列，令S的其余列为零，得到稀疏矩阵S。

示例性的，通过以上子步骤3c)-3h)，迭代求解稀疏矩阵S的过程如下所示：

第一次迭代，i＝1：

I)设定初值P₀＝Z，利用公式 $P_0=Z-{\mathrm{\Φ}}_{{\mathrm{\Ω}}_0}{\left({{\mathrm{\Φ}}^H}_{{\mathrm{\Ω}}_0}{\mathrm{\Φ}}_{{\mathrm{\Ω}}_0}\right)}^{-1}{{\mathrm{\Φ}}^H}_{{\mathrm{\Ω}}_0}Z$ 求出其中Ω₀＝Φ；

II)利用公式 ${\mathrm{\β}}_1=\arg \min \left[\mathrm{\Σ}\right|({\mathrm{\Φ}}_{{\mathrm{\Ω}}_0}^H{P}_0-{\left|\right|{\mathrm{\Φ}}_{{\mathrm{\Ω}}_0}^H{P}_0\left|\right|}_{\∞})\left|\right],$ 求出β₁；

III)利用公式Ω₁＝Ω₀∪β₁求出Ω₁；

第二次迭代，i＝2：

I)利用第一次迭代求出的Ω₁，求出利用公式 $P_1=Z-{\mathrm{\Φ}}_{{\mathrm{\Ω}}_1}{\left({{\mathrm{\Φ}}^H}_{{\mathrm{\Ω}}_1}{\mathrm{\Φ}}_{{\mathrm{\Ω}}_1}\right)}^{-1}{{\mathrm{\Φ}}^H}_{{\mathrm{\Ω}}_1}Z$ 求出P₁；

II)利用公式 ${\mathrm{\β}}_2=\arg \min \left[\mathrm{\Σ}\right|({\mathrm{\Φ}}_{{\mathrm{\Ω}}_1}^H{P}_1-{\left|\right|{\mathrm{\Φ}}_{{\mathrm{\Ω}}_1}^H{P}_1\left|\right|}_{\∞})\left|\right],$ 求出β₂；

III)利用公式Ω₂＝Ω₁∪β₂求出Ω₂

同理第K步迭代得到Ω_K。设定Ω＝Ω_K，Ω_K＝Ω_K-1∪β_k共有K个元素。

3k)将阵列流型矩阵D和稀疏矩阵S代入优化模型为(10)中，从而估计出N个阵元的扰动矩阵Γ₁(θ)：

${\mathrm{\Γ}}_1\left(\mathrm{\θ}\right)={\left|\right|Z-[A\left(\mathrm{\θ}\right)+D]S\left|\right|}_F^2+\mathrm{\μ}{\left|\right|S\left|\right|}_{\∞,0}---\left(14\right)$

其中，其中，为Frobenius范数运算符，Z表示回波信号的模型矩阵，A(θ)是目标的导向矢量，S表示稀疏矩阵，D表示阵列流行矩阵；||·||_∞，0表示混合范数，μ>0为正则化参数，θ表示目标的到达角。

现有技术求解阵元的扰动矩阵Γ₁(θ)时都会有一定的偏差；而本发明在步骤3中通过建立最优化模型然后分别求出阵列流行矩阵D和稀疏矩阵S，进而求出阵元的扰动矩阵Γ₁(θ)，以便利用Γ₁(θ)估计出阵元的位置偏差矩阵(ΔX₁,ΔY₁)，以便估计出阵元的位置(X₁,Y₁)。

步骤4，利用估计出的N个阵元的扰动矩阵Γ₁(θ)，估计出N个阵元的位置偏差矩阵(ΔX₁,ΔY₁)，(ΔX₁,ΔY₁)＝[(ΔX₁,ΔY₁)₁,(ΔX₁,ΔY₁)₁,...,((ΔX₁,ΔY₁)_N)]^T。

将步骤3求出的阵元的扰动矩阵Γ₁(θ)代入(15)式中，求出N个阵元的位置偏差矩阵(X₁,Y₁)：

$\left\{\begin{array}{c}{P}_1\left(\mathrm{\θ}\right)=-\frac{\mathrm{\λ}}{2\mathrm{\π}}\mathrm{angle}\left({\mathrm{\Γ}}_1\left(\mathrm{\θ}\right)\right)\\ ({\mathrm{\ΔX}}_1,{\mathrm{\ΔY}}_1)=\left[\begin{array}{cc}{P}_1\left({\mathrm{\θ}}_1\right)& P_1\left({\mathrm{\θ}}_2\right),...,P_1\left({\mathrm{\θ}}_l\right),...,P_1\left({\mathrm{\θ}}_K\right)\end{array}\right]{\left({\left[\begin{array}{cccccc}{\sin \mathrm{\θ}}_1& {\sin \mathrm{\θ}}_2& & {\sin \mathrm{\θ}}_l& & {\sin \mathrm{\θ}}_K\\ & & ...& & ...& \\ {\cos \mathrm{\θ}}_1& {\cos \mathrm{\θ}}_2& & {\cos \mathrm{\θ}}_l& & {\cos \mathrm{\θ}}_K\end{array}\right]}^T\right)}^{-1}\end{array}\right.---\left(15\right)$

其中，λ表示雷达的波长，[]^-1表示求矩阵的逆，Γ₁(θ)表示阵元的扰动矩阵，angle为求向量的角度，θ表示目标的到达角。

步骤5，利用估计出的N个阵元的位置偏差矩阵(ΔX₁,ΔY₁)，得到估计出的阵元位置(X₁,Y₁)，估计出的阵元位置(X₁,Y₁)＝(X',Y')+(ΔX₁,ΔY₁)。

下面结合仿真实验对本发明的效果做进一步说明。

仿真条件：本发明的仿真是在MATLAB R2009a的软件环境下进行的。阵元的幅度扰动矩阵A(θ)的均值为1、方差为0.2，阵元的相位扰动矩阵的均值为0、方差为2，干扰矩阵E均值为0，方差为1。雷达天线的阵元数为N＝13。目标的到达角θ为：θ₁＝10°和θ₂＝30°，雷达的波长λ设为0.3米。

本发明中雷达天线设为均匀线阵，且倾斜放置。10个阵元的位置矩阵(X',Y')分别为：(2,3)、(2.1,3.15)、(2.2，3.3)、(2.3,3.45)、(2.4,3.6)、(2.5,3.75)、(3.6,3.9)、(2.7,4.05)、(2.8,4.2)和(2.9,4.35)，即，阵元间的直线距离为0.15(λ/2)米。

仿真内容1：用DBF、MUSIC和自校正算法估计到达角θ。

仿真结果：如图2(a)和2(b)所示，实线表示自校正算法计算出的到达角θ，虚线表示DBF算法计算出的到达角θ，点画线表示MUSIC算法计算出的到达角θ。横坐标表示到达角θ单位为度，纵坐标为幅度，单位为dB。可以发现用DBF和MUSIC估计出的到达角θ与设定的到达角θ会有一定的偏差，而用自校正算法估计出的到达角θ与设定的到达角θ相一致。所以自校正算法能比较精确地估计出到达角θ。本发明方案中步骤2采用自校正算法估计目标的到达角。由仿真图1可以得到：用DBF得到到达角θ为：θ₁＝8°，θ₂＝28°；用MUSIC得到的到达角θ为：θ₁＝8°，θ₂＝28°；用自校正算法得到的到达角θ为：θ₁＝10°，θ₂＝30°，与设定的θ相同。

仿真内容2：将初始扰动矩阵Γ(θ)代入(15)中，计算出实际的阵元的位置偏差矩阵(ΔX,ΔY)并和用本专利提出的方法估计出的阵元位置偏差矩阵(ΔX₁,ΔY₁)进行比较分析。

仿真结果：如图3所示，圈形表示实际的阵元位置偏差(ΔX,ΔY)，星形表示用本发明估计出的阵元位置偏差(ΔX₁,ΔY₁)。横坐标表示阵元数，纵坐标表示阵元的位置偏差，单位为米，可以发现实际的阵元位置偏差和用本专利提出的方法估计的阵元的位置偏差很接近。所以，本专利提出的方法，可以比较准确的估计出阵元的位置偏差，有利于对阵元的位置进行校正。

仿真内容3：用实际的位置偏差矩阵(ΔX,ΔY)和用本发明提出的方法估计出的阵元位置偏差矩阵(ΔX₁,ΔY₁)分别对阵元的位置矩阵(X',Y')进行校正。

仿真结果：如图4所示，横坐标表示阵元的水平位置，纵坐标表示阵元的高度。

用实际的阵元的位置偏差矩阵(ΔX,ΔY)与阵元的位置矩阵(X',Y')相加，得到阵元的实际位置矩阵(X,Y)为：(2.1426,3.1426)、(2.1537,3.2037)、(2.2531,3.3531)、(2.4530,3.6030)、(2.4852,3.6852)、(2.6154,3.8654)、(2.6816,3.9816)、(2.7585,4.1085)、(2.8142,4.2142)和(2.9150,4.3650)，图4中，圈形表示阵元的实际位置。

本发明估计出的位置偏差矩阵(ΔX₂,ΔY₂)与阵元的位置矩阵(X',Y')相加，得到估计出的阵元的位置矩阵(X₂,Y₂)为：(2.1420,3.1420)、(2.1540,3.2040)、(2.2522,3.3522)、(2.4529,3.6029)、(2.4847,3.6847)、(2.6142,3.8642)、(2.6833,3.9833)、(2.7574,4.1074)、(2.8135,4.2135)和(2.9137,4.3637)，图4中星形表示估计出的阵元的位置。

由图4可以看出，阵元的实际位置和用本发明估计出的阵元位置基本相同，所以本发明能够比较精确地估计出阵元的位置。

Claims (5)

1.一种雷达天线阵元位置的校正方法，其特征在于，包括以下步骤：

步骤1，雷达天线接收原始的回波信号，建立该回波信号的模型矩阵Z；利用回波信号的模型矩阵Z求出回波信号的自相关矩阵R；并且设定雷达天线的N个阵元的位置矩阵(X',Y')，其中(X',Y')＝[(X'₁,Y'₁)(X'₂,Y'₂)...(X'_N,Y'_N)]^T表示N个阵元的位置矩阵；

步骤2，利用回波信号的模型矩阵Z和回波信号的自相关矩阵R估计出目标的到达角θ；构造到达角相关矩阵Q(θ)，并求出到达角相关矩阵Q(θ)的特征值对应的特征向量构成矩阵V；

步骤3，根据目标的到达角θ和到达角相关矩阵Q(θ)的特征值对应的特征向量构成矩阵V，估计出N个阵元的扰动矩阵Γ₁(θ)，Γ₁(θ)＝[(Γ₁(θ))₁,(Γ₁(θ))₂,...,(Γ₁(θ))_N)]^T；

步骤4，利用估计出的N个阵元的扰动矩阵Γ₁(θ)，估计出N个阵元的位置偏差矩阵(ΔX₁,ΔY₁)，(ΔX₁,ΔY₁)＝[(ΔX₁,ΔY₁)₁,(ΔX₁,ΔY₁)₁,...,((ΔX₁,ΔY₁)_N)]^T；

步骤5，利用估计出的N个阵元的位置偏差矩阵(ΔX₁,ΔY₁)，得到估计出的阵元位置(X₁,Y₁)，估计出的阵元位置(X₁,Y₁)＝(X',Y')+(ΔX₁,ΔY₁)。

2.根据权利要求1所述的一种雷达天线阵元位置的校正方法，其特征在于，步骤1包括以下子步骤：

1a)构建目标的导向矢量矩阵A(θ)：

$A\left(\mathrm{\θ}\right)=\exp (-j\frac{2\mathrm{\π}}{\mathrm{\λ}}\×(0:N-1)\×d\×\sin (\mathrm{\θ}\×\frac{\mathrm{\π}}{180}))---\left(1\right)$

其中,N表示雷达天线的阵元数，λ表示雷达工作的波长，d表示阵元的间隔距离,θ表示目标的到达角，θ＝[θ₁,θ₂,...θ_l,...,θ_K],θ_l表示第l个目标的到达角，K表示目标的个数；

1b)利用目标的导向矢量矩阵A(θ)，构建雷达天线的N个阵元的初始扰动矩阵Γ(θ)：

其中，Γ(θ)为N×1维矩阵，N为雷达天线的阵元数，A(θ)表示目标的导向矢量矩阵，表示阵元相位的随机扰动矩阵，θ表示目标的到达角；

1c)利用目标的导向矢量矩阵A(θ)和初始扰动矩阵Γ(θ)构造回波信号的模型矩阵Z；

Z＝Γ(θ)A(θ)+E      (3)

其中，A(θ)表示目标的导向矢量矩阵，Γ(θ)表示N个阵元的初始扰动矩阵，E表示噪声矩阵，θ表示目标的到达角；

利用回波信号的模型矩阵Z求出回波信号的自相关矩阵R：

$R=\frac{1}{N}{\mathrm{ZZ}}^H---\left(4\right)$

其中，Z表示回波信号的模型矩阵，N表示雷达天线的阵元数，[·]^H表示共轭转置；

1d)设定雷达天线的N个阵元的位置矩阵(X',Y')，其中，(X',Y')＝[(X₁',Y₁'),(X₂',Y₂'),...,(X_N',Y_N')]^T。

3.根据权利要求1所述的一种雷达天线阵元位置的校正方法，其特征在于，步骤2包括以下子步骤：

2a)对目标的自相关矩阵R进行特征值分解，得到目标子空间U_s和噪声子空间U_N；

2b)利用目标的导向矢量矩阵A(θ)和阵元的初始扰动矩阵Γ(θ)构造目标辅助矩阵w(θ)；

w(θ)＝a_m(θ)Γ(θ)      (5)

其中，a_m(θ)由A(θ)的第m列构成，m＝1,2,3...,K，K表示目标的个数，θ表示目标的到达角；

2c)根据目标辅助矩阵w(θ)与噪声子空间U_N的正交性，得到下式：

w^H(θ)U_NU^H _Nw(θ)＝0      (6)

其中，U_N表示噪声子空间，N表示雷达天线阵元数，[·]^H表示共轭转置，θ表示目标的到达角；

将目标辅助矩阵w(θ)公式(5)代入公式(6)得：

Γ^H(θ)a_m ^H(θ)U_NU^H _Na_m(θ)Γ(θ)＝0      (7)

由公式(7)得到达角相关矩阵Q(θ)：

Q(θ)＝a_m ^H(θ)U_NU^H _Na_m(θ)      (8)

2d)将到达角相关矩阵Q(θ)代入以下公式(9)求出目标的到达角θ：

$\mathrm{\θ}=\arg \underset{\mathrm{\θ}}{\max }\frac{1}{{\mathrm{\λ}}_{\min }\left[Q\left(\mathrm{\θ}\right)\right]}$ 或 $\mathrm{\θ}=\arg \underset{\mathrm{\θ}}{\max }\frac{1}{\det \left[Q\left(\mathrm{\θ}\right)\right]}---\left(9\right)$

其中，arg[·]为求解最优化，max表示求最大值，λ_min[·]为求矩阵的最小特征值，[·]^H表示矩阵的共轭转置，det[·]为求矩阵的行列式，θ表示目标的到达角；

2e)对到达角相关矩阵Q(θ)进行奇异值分解，得到达角相关矩阵Q(θ)的特征值对应的特征向量构成矩阵V。

4.根据权利要求1所述的一种雷达天线阵元位置的校正方法，其特征在于，步骤3包括以下子步骤：

3a)将求解阵元的扰动矩阵Γ₁(θ)转换为如下的优化模型公式(10)：

${\mathrm{\Γ}}_1\left(\mathrm{\θ}\right)=\underset{S.D}{\min }{\left|\right|Z-[A\left(\mathrm{\θ}\right)+D]S\left|\right|}_F^2+\mathrm{\μ}{\left|\right|S\left|\right|}_{\∞,0}---\left(10\right)$

其中，||·||_F为Frobenius范数运算符，Z表示回波信号的模型矩阵，A(θ)是目标的导向矢量，S为目标的稀疏矩阵，D为阵列流行矩阵；||·||_∞，0表示混合范数，μ>0为正则化参数，θ表示目标的到达角；

3b)设定阵列流型矩阵D等于Q(θ)的特征值对应的特征向量构成矩阵V，即：

D＝V      (11)

3c)将优化模型公式(10)转换为以下模型公式(12)：

$\underset{S}{\min }{\left|\right|Z-[A\left(\mathrm{\θ}\right)+D]S\left|\right|}_F^2---\left(12\right)$

其中，Z表示回波信号的模型矩阵，A(θ)表示目标的导向矢量矩阵，S表示Q(θ)的非零行构成的稀疏矩阵，D表示阵列流行矩阵，D等于到达角相关矩阵Q(θ)的特征值对应的特征向量构成矩阵V，θ表示目标的到达角；

3d)设定Φ＝A(θ)D，然后设定Φ的支撑矩阵Ω，支撑矩阵Ω中包含了Φ中不为零的列；

3e)β是支撑矩阵Ω中的一个元素，通过求解式(13)得到β的第i次迭代元素β_i：

${\mathrm{\β}}_i=\arg \min \left[\mathrm{\Σ}\right|({\mathrm{\Φ}}_{{\mathrm{\Ω}}_{i-1}}^H{P}_{i-1}-{\left|\right|{\mathrm{\Φ}}_{{\mathrm{\Ω}}_{i-1}}^H{P}_{i-1}\left|\right|}_{\∞})\left|\right]---\left(13\right)$

其中，表示在第i-1次迭代中由支撑矩阵Ω_i-1中的元素对应于矩阵Φ中的列向量构成的矩阵，设定支撑矩阵Ω的初值Ω₀＝Φ；

构造P_i-1表示矩阵P第i-1次迭代矩阵，设定矩阵P的初值P₀＝Z，[·]^H表示共轭转置，arg表示求解最优化问题，i＝1,2,...,K，K表示目标数；

3f)利用第i次迭代元素β_i求解第i次迭代时的支撑矩阵Ω_i＝Ω_i-1∪β_i；

3g)令迭代次数i增加1，重复迭代3e)-3f)步，直到i等于K，得到Ω_K迭代结束；K表示目标的个数，并设定支撑矩阵Ω＝Ω_K，Ω_K＝Ω_K-1∪β_k共有K个元素；

3h)利用求出的支撑矩阵Ω，求出稀疏矩阵S：

定义矩阵S_Ω，S_Ω计算公式为其中，Φ_Ω表示由支撑矩阵Ω中的元素对应于矩阵Φ中的列向量构成的矩阵，[·]^H表示共轭转置操作，Z表示回波信号的模型矩阵，[·]^-1表示求矩阵的逆，A(θ)是目标的导向矢量，θ表示目标的到达角；

矩阵S_Ω组成目标的稀疏矩阵S的非零列，令S的其余列为零，得到稀疏矩阵S；

3k)将阵列流型矩阵D和稀疏矩阵S代入优化模型为(10)中，估计出N个阵元的扰动矩阵Γ₁(θ)：

${\mathrm{\Γ}}_1\left(\mathrm{\θ}\right)={\left|\right|Z-[A\left(\mathrm{\θ}\right)+D]S\left|\right|}_F^2+\mathrm{\μ}{\left|\right|S\left|\right|}_{\∞,0}---\left(14\right)$

其中，其中，为Frobenius范数运算符，Z表示回波信号的模型矩阵，A(θ)是目标的导向矢量，S表示稀疏矩阵，D表示阵列流行矩阵，D等于到达角相关矩阵Q(θ)的特征值对应的特征向量构成矩阵V；||·||_∞，0表示混合范数，μ>0为正则化参数，θ表示目标的到达角。

5.根据权利要求1所述的一种雷达天线阵元位置的校正方法，其特征在于，步骤4具体包括：

将步骤3估计的阵元的扰动矩阵Γ₁(θ)代入以下公式(15)中，求出N个阵元的位置偏差矩阵(ΔX₁,ΔY₁)：

$\left\{\begin{array}{c}{P}_1\left(\mathrm{\θ}\right)=-\frac{\mathrm{\λ}}{2\mathrm{\π}}\mathrm{angle}\left({\mathrm{\Γ}}_1\left(\mathrm{\θ}\right)\right)\\ ({\mathrm{\ΔX}}_1,{\mathrm{\ΔY}}_1)=\left[\begin{array}{cc}{P}_1\left({\mathrm{\θ}}_1\right)& P_1\left({\mathrm{\θ}}_2\right),...,P_1\left({\mathrm{\θ}}_l\right),...,P_1\left({\mathrm{\θ}}_K\right)\end{array}\right]{\left({\left[\begin{array}{cccccc}{\sin \mathrm{\θ}}_1& {\sin \mathrm{\θ}}_2& & {\sin \mathrm{\θ}}_l& & {\sin \mathrm{\θ}}_K\\ & & ...& & ...& \\ {\cos \mathrm{\θ}}_1& {\cos \mathrm{\θ}}_2& & {\cos \mathrm{\θ}}_l& & {\cos \mathrm{\θ}}_K\end{array}\right]}^T\right)}^{-1}\end{array}\right.---\left(15\right)$

其中，λ表示雷达的波长，[]^-1表示求矩阵的逆，Γ₁(θ)表示估计出的阵元的扰动矩阵，angle为求向量的角度，θ表示目标的到达角，K表示目标的个数。