一种稀布阵列近场无源定位幅相误差校正方法

摘要

该发明公开了一种稀布阵列近场无源定位幅相误差校正方法,涉及相控阵雷达技术领域。在目标位于近场的情况下，固定俯仰角将目标距离和方位角进行交替循环估计，直至两个参数都收敛至真实值，作为误差校正的初始位置；然后根据，对目标方位角进行一维搜索并用上述幅相误差自校正的方法对阵列幅相误差进行校正，估计出阵列幅相误差矩阵Γθ以及方位角的一次迭代值θ′；估计出阵列幅相误差矩阵Γr以及目标距离的一次迭代值；以此进行循环迭代，直至||Гθ-Γr||2，即参数估计值收敛于真实值。本发明在一定的布阵范围内，阵元位置随机稀疏摆放，只需要较少的阵元就可达到相同的角度分辨率，在实际工程中有利于节约成本。

说明书

技术领域

本发明涉及相控阵雷达技术领域，具体涉及无源定位的阵列误差校正方法。

背景技术

与常规均匀线阵相比，稀布阵列用相同的阵元使得阵列孔径增大，以较少的阵元 达到较高的分辨率，从而简化阵列结构，降低造价。目前，稀布阵列天线在抗环境干扰的卫 星接收天线，高频地面雷达天线和射电天文中的干涉阵列等领域中正在得到越来越广泛的 应用。另外，波达方向(DirectionofArrival,DOA)估计在雷达等领域广泛应用，通常高分 辨DOA估计阵列处理均基于理想阵列信号模型，但是具体在实际工程应用中，各通道间硬件 参数的不同，使得目标回波信号在不同的阵元位置有不同的幅度和相位加权，即阵列中存 在多通道幅相误差问题。

WeissA.J.和FriedlanderB.等人在“Eigenstructuremethodsfordirection findingwithsensorgainandphaseuncertainties”(Acoustics,Speech,andSignal Processing,1988.ICASSP-88.,1988InternationalConferenceon.IEEE,1988:2681- 2684)中提出了将阵列通道幅相误差参数与信号源方位交替联合迭代的阵列误差自校正算 法(称为WF算法)，以上文献只针对了在远场情况下对于MUSIC算法进行的幅相误差估计，但 是在信源位于近场的情况下，信号的波前不能近似为平面波而要看作球面波，实现定位需 要在距离和角度上进行联合二维搜索。本发明提供了一种基于迭代的稀布阵列近场无源定 位幅相误差校正方法，在距离维和角度维上进行轮换迭代，即使在幅相误差较大时，该方法 也能精确校正稀布阵列通道幅相误差。

发明内容

为了克服上述现有技术的不足，本发明提供了一种稀布阵列近场无源定位幅相误 差校正方法,从而达到定位精度高，成本低廉的目的。

本发明提出的基于迭代的稀布阵列近场无源定位福相误差校正方法基本思想如 下：在目标位于近场的情况下，固定俯仰角将目标距离r和方位角θ进行交替循环估计， 直至两个参数都收敛至真实值，即：首先在存在误差的情况下粗略估计出目标距离r₀和方 位角θ₀，作为误差校正的初始位置；然后根据r₀，对目标方位角θ进行一维搜索并用上述幅相 误差自校正的方法对阵列幅相误差进行校正，估计出阵列幅相误差矩阵Γ_θ以及方位角θ的 一次迭代值θ′；最后根据θ′，对目标距离r进行一维搜索并用上述幅相误差自校正的方法对 阵列幅相误差进行校正，估计出阵列幅相误差矩阵Γ_r以及目标距离r的一次迭代值r′；以 此进行循环迭代，直至||Γ_θ-Γ_r||<ε₂，即参数估计值收敛于真实值，迭代过程如附图2所 示。因而本发明是一种稀布阵列近场无源定位幅相误差校正方法，该方法具体包括：

步骤1：获得定位场中信号接收机即定位阵元的个数及其位置信息；

步骤2：以第一个阵元为参考阵元，各阵元接收P个位于定位场中的目标信号源发 送的信号；

步骤3：用x_m(t)表示第m个阵元在第t个时刻接收信号的采样值，则天线阵列接收 信号

$\left(\begin{array}{cc}x\left(t\right)={[x_1\left(t\right),x_2\left(t\right),...,x_M\left(t\right)]}^T=\underset{p=1}{\overset{P}{\Sigma }}a(r_p,{\theta }_p)s_p\left(t\right)+n\left(t\right),& t=1,2,...,L\end{array}\right)$

其中，(·)^T表示转置操作符，s_p(t)表示第p个近场信号源，n(t)表示天线阵列接收 的噪声矢量，L表示采样总数，a(r_p,θ_p)表示第p个近场信号源导向矢量，r_p表示第p个近场信 号源到参考阵元的距离；

步骤4：通过接收天线阵列数据x(t)估算协方差矩阵(·)^H表 示对矩阵进行共轭转置计算；再估算协方差矩阵的逆矩阵R^-1，(·)^-1表示对矩阵进行求逆 运算；

步骤5：采用最小方差波束形成算法估计出目标角度和距离的初始值 ${\hat{\theta }}^{\left(0\right)}=({\hat{\theta }}_1^{\left(0\right)},{\hat{\theta }}_2^{\left(0\right)},...,{\hat{\theta }}_P^{\left(0\right)}),$令二维迭代次数N＝0；

步骤6：首先在角度维上进行迭代，固定目标距离为则可以在角度维上进行一 维幅相误差校正，设定角度维上迭代次数变量为k，从角度维上估计出的稀布阵列雷达幅相 误差对角矩阵为初始值为

定义代价函数

$J_k=\underset{p=1}{\overset{P}{\Sigma }}{a}^H({r_p}^{\left(k\right)},{{\theta }_p}^{\left(k\right)}){\left({{\Gamma }_{\theta }}^{\left(k\right)}\right)}^H{R}^{-1}{{\Gamma }_{\theta }}^{\left(k\right)}a({r_p}^{\left(k\right)},{{\theta }_p}^{\left(k\right)})$

其中表示第p个近场信号源的导向矢量，(·)^H表示对矩阵进行共轭 转置计算；增加迭代次数，直到J_k-1-J_k≤ε₁结束迭代，获得角度估计值和幅相误差估计 值为

步骤7：在距离维上进行迭代，固定目标角度为则可以在距离维上进行一维 幅相误差校正；设定从距离维上估计出的稀布阵列雷达幅相误差对角矩阵初始值为 进行与步骤7相同的一维迭代过程，可求得距离估计值和幅相误差估计值

步骤8：比较两个维度上估计出的幅相误差对角矩阵，若满足则停止 迭代，其中ε₂为根据实际情况预设的门限，即得到幅相误差的估计值否则，令N＝N+ 1，转到步骤6，继续迭代，其中N为角度维和距离维共同迭代次数的变量。

进一步的，所述步骤6的具体步骤为：

步骤6.1：定义空间谱公式为:

$p\left({\hat{\theta }}^{\left(k\right)}|{{\hat{\Gamma }}_{\theta }}^{\left(k\right)},{\hat{r}}^{\left(k\right)}\right)=\frac{1}{{\hat{a}}^H\left({\hat{r}}^{\left(k\right)},{\hat{\theta }}^{\left(k\right)}\right){\left({{\hat{\Gamma }}_{\theta }}^{\left(k\right)}\right)}^H{R}^{-1}{{\hat{\Gamma }}_{\theta }}^{\left(k\right)}\hat{a}\left({\hat{r}}^{\left(k\right)},{\hat{\theta }}^{\left(k\right)}\right)}$

根据利用上述空间谱公式求取对应的P个最大值时的角度便可得到第k次迭代的P个波达方向的估计值；

步骤6.2：令w＝[1,0,...,0]^T，则满足约束条件则用拉格朗日乘子 法，可求得幅相误差矩阵如下:

$\left(\begin{array}{c}{{\hat{\Gamma }}_{\theta }}^{(k+1)}=diag\left\{{\hat{\delta }}^{(k+1)}\right\}\\ {\hat{\delta }}^{(k+1)}=\frac{{\left(Q_k\right)}^{-1}w}{w^T{\left(Q_k\right)}^{-1}w}\end{array}\right)$

其中，$Q_k=\underset{p=1}{\overset{P}{\Sigma }}{\tilde{a}}^H({\hat{r}}_p^{\left(k\right)},{\hat{\theta }}_p^{\left(k\right)})R^{-1}\tilde{a}(r_p^{\left(k\right)},{\hat{\theta }}_p^{\left(k\right)}),\tilde{a}({\tilde{r}}_p^{\left(k\right)},{\hat{\theta }}_p^{\left(k\right)})=diag\left\{a({\hat{r}}_p^{\left(k\right)},{\hat{\theta }}_p^{\left(k\right)})\right\};$

步骤6.3判断是否收敛

$J_k={\left({\hat{\delta }}^{(k+1)}\right)}^H{Q}_k{\hat{\delta }}^{(k+1)}$

在满足J_k-1-J_k>ε₁，时,继续迭代，ε₁为根据实际情况预设的门限值；否则,结束迭 代；最后一次迭代得到角度估计值和幅相误差估计值为

与现有技术相比，本发明方法能够对近场目标进行阵列幅相误差的二维校正，并 利用近场无源定位的二维特点，对估计出的误差矩阵进行修正并能实现准确估计。此外，本 发明基于一种随机稀布阵列，即在一定的布阵范围内，阵元位置随机稀疏摆放，只需要较少 的阵元就可达到相同的角度分辨率，在实际工程中有利于节约成本。

附图说明

图1为近场窄带信号无源定位模型；

图2为一种基于迭代的稀布阵列近场无源定位幅相误差校正方法流程图；

图3为存在阵列误差情况下，稀布阵列近场无源定位的方向图；

图4为存在阵列误差情况下，采用本发明的幅相误差校正方法，对阵列接收信号进 行补偿后，稀布阵列近场无源定位的方向图；

图5是图3和图4中当方位角为-22.5°时，误差校正前后的切面对比图；

图6是图3和图4中当距离为1700m时，误差校正前后的切面对比图。

具体实施方式

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

步骤1:在笛卡尔坐标系中的x轴上的[0,D]米范围内随机布置M个天线，形成一个 稀布阵，其位置产生方式如下：在位置0以及位置D处各布置1个阵元且位置分别设为d₁与d_M， 即d₁＝0，d_M＝D，固定阵列孔径。然后在(0,D)范围以均匀分布方式随机产生M-2个数，并从小 到大进行排序，将剩下的M-2个阵元分别布置在M-2个随机数表示的位置上，并设阵元位置 分别为d₂,…,d_m,…，d_M-1。

步骤2:以d₁所在的位置的阵元为参考阵元，考虑P个近场不相关目标信号源入射 到步骤1产生的阵列，对任意的p＝1,2,…,P，用(r_p,θ_p)描述第p个信号源的位置，其中，如图 1所示，r_p表示第p个信号源到参考阵元的距离，θ_p表示第p个信号源相对于参考阵元与z轴的 夹角；此外，P个近场不相关目标信号源要满足如下条件；

近场条件：0<r_p≤2D²/λ_p，对任意的p＝1,2,…,P，其中，λ_p表示第p个信号源的波 长；

不相关条件：

$\left(\begin{array}{c}r==\frac{E\left\{\right[s_i-E\left\{s_i\right\}\left]\right[s_q-E\left\{s_q\right\}\left]\right\}}{\sqrt{E\left\{{[s_i-E\{s_i\left\}\right]}^2\right\}E\left\{{[s_q-E\{s_q\left\}\right]}^2\right\}}}=0\\ \begin{array}{ccc}i=1,2,\mathrm{...},P& q=1,2,\mathrm{...},P& i\ne q\end{array}\end{array}\right)$

其中，E[·]表示求期望，s_i表示第i个近场源信号，s_q表示第q个近场源信号，r表示 相关系数，对任意的i＝1,2,…,Pq＝1,2,…,P且i≠q，r＝0表示s_i与s_q不相关。

步骤3：则由几何知识可得第p个信号源到第m个阵元的距离r_mp可表示为

$r_{mp}=\sqrt{{(r_p{\mathrm{sin\theta }}_p-d_m)}^2+{\left(r_p{\mathrm{cos\theta }}_p\right)}^2}$

则第p个信号源到第m个阵元与第p个信号源到参考阵元的距离差△r_mp可以表示为 如下公式

$\left(\begin{array}{c}{\mathrm{\Delta r}}_{mp}=\sqrt{{\left(r_p{\mathrm{sin\theta }}_p-d_m\right)}^2+{\left(r_p{\mathrm{cos\theta }}_p\right)}^2}-r_p\\ =r_p{\left(1+\frac{d_m^2}{r_p^2}-2\frac{d_m{\mathrm{sin\theta }}_p}{r_p}\right)}^{\frac{1}{2}}-r_p\end{array}\right)$

若当距离参数满足0.62(D³/λ_p)^1/2<r_p<2D²/λ_p，即第p个近场信源位于菲涅耳区，则 距离差△r_mp可由泰勒级数展开近似表示为

${\mathrm{\Delta r}}_{mp}=-d_m{\mathrm{sin\theta }}_p+\frac{d_m^2{\cos }^2{\theta }_p}{2r_p}$

则第m个阵元在第t个采样时刻接收的数据表示为

$\left(\begin{array}{cc}{x}_m\left(t\right)=\underset{p=1}{\overset{P}{\Sigma }}\frac{r_p}{r_{mp}}{s}_p\left(t\right)e^{-j\frac{2\pi }{{\lambda }_p}\left(-d_m{\mathrm{sin\theta }}_p+\frac{d_m^2{\cos }^2{\theta }_p}{2r_p}\right)}+n\left(t\right),& t=1,2,...,L\end{array}\right)$

其中，r_mp表示第p个信号源到第m个阵元的距离，n(t)表示t时刻噪声信号，L表示总 的采样数。

将M个阵元在t时刻采样的数据排成一个列矢量x(t)，可以表示为

$x\left(t\right)={[x_1\left(t\right),x_2\left(t\right),...,x_M\left(t\right)]}^T=\underset{p=1}{\overset{P}{\Sigma }}a(r_p,{\theta }_p)s_p\left(t\right)+n\left(t\right),$

其中，n(t)是独立的高斯白噪声信号矢量，a(r_p,θ_p)表示第p个近场信号源导向矢 量，其具体表达式为

$a\left(r_p,{\theta }_p\right)={\left[1,\frac{r_p}{r_{2p}}{e}^{-j\frac{2\pi }{{\lambda }_p}\left(-d_2{\mathrm{sin\theta }}_p+\frac{d_2^2{\cos }^2{\theta }_p}{2r_p}\right)},\mathrm{...},\frac{r_p}{r_{Mp}}{e}^{-j\frac{2\pi }{{\lambda }_p}\left(-d_M{\mathrm{sin\theta }}_p+\frac{d_M^2{\cos }^2{\theta }_p}{2r_p}\right)}\right]}^T$

步骤4：估计接收数据的协方差矩阵可以表示为

$\hat{R}=\frac{1}{L}\underset{t=1}{\overset{L}{\Sigma }}x\left(t\right)x^H\left(t\right)$

计算协方差矩阵的逆其中(·)^-1表示对矩阵求逆。

步骤5:设定迭代次数变量为k，第k次迭代中从距离维上估计出的稀布阵列雷达幅 相误差对角矩阵为从角度维上估计出的稀布阵列雷达幅相误差对角矩阵为令初 始值k＝0，取其中Γ₀表示雷达幅相误差对角矩阵初始值。

步骤6:进行迭代之前先结合最小方差波束形成算法估计出目标角度和距离的初 始值

${\hat{r}}^{\left(0\right)}=({\hat{r}}_1^{\left(0\right)},{\hat{r}}_2^{\left(0\right)},...,{\hat{r}}_p^{\left(0\right)}),{\hat{\theta }}^{\left(0\right)}=({\hat{\theta }}_1^{\left(0\right)},{\hat{\theta }}_2^{\left(0\right)},\mathrm{...},{\hat{\theta }}_P^{\left(0\right)}).$

最小方差波束形成算法的基本思想是：

通过对接收信号进行加权求和使得期望方向的功率最强，从而实现目标方向信号 最大化接收并抑制其他方向干扰，其期望方向的加权矢量跟方向矢量相同。最小均方误差 准则可以总结

为：

从数学上可以理解，最小方差算法就是在满足的所以权向量w条 件下，找到满足使得阵列输出的功率w^HR_xw最小的权向量w。从物理意义上理解，上式保证了 θ₀方向上的信号正确接收，抑制了其它方向上的干扰信号。

对上式采用拉格朗日乘子法可得最优的权矢量为

则阵列的谱函数为

对谱函数进行搜索，得到的谱峰值就是来波方向以及位置。

步骤7:首先在角度维上进行迭代，固定目标距离为则可以在角度维上进行一 维幅相误差校正。

定义代价函数

$J_k=\underset{p=1}{\overset{P}{\Sigma }}{a}^H({{\theta }_p}^{\left(k\right)},{\hat{r}}^{\left(k\right)}){\left({{\Gamma }_{\theta }}^{\left(k\right)}\right)}^H{R}^{-1}{{\Gamma }_{\theta }}^{\left(k\right)}a({{\theta }_p}^{\left(k\right)},{\hat{r}}^{\left(k\right)})$

式中：“||·||”表示Frobenius范数。

7.1定义空间谱公式为:

$p\left({\hat{\theta }}^{\left(k\right)}|{{\hat{\Gamma }}_{\theta }}^{\left(k\right)},{\hat{r}}^{\left(k\right)}\right)=\frac{1}{{\hat{a}}^H\left({\hat{\theta }}^{\left(k\right)},{\hat{r}}^{\left(k\right)}\right){\left({{\hat{\Gamma }}_{\theta }}^{\left(k\right)}\right)}^H{R}^{-1}{{\hat{\Gamma }}_{\theta }}^{\left(k\right)}\hat{a}\left({\hat{\theta }}^{\left(k\right)},{\hat{r}}^{\left(k\right)}\right)}$

根据利用上述空间谱公式求取对应的P个最大值时的角度便 可得到第k次迭代的P个波达方向的估计值。

7.2令w＝[1,0,...,0]^T，则满足约束条件则利用拉格朗日乘子法， 可求得幅相误差矩阵如下:

$\left(\begin{array}{c}{{\hat{\Gamma }}_{\theta }}^{(k+1)}=diag\left\{{\hat{\delta }}^{(k+1)}\right\}\\ {\hat{\delta }}^{(k+1)}=\frac{{\left(Q_k\right)}^{-1}w}{w^T{\left(Q_k\right)}^{-1}w}\end{array}\right)$

其中，$Q_k=\underset{p=1}{\overset{P}{\Sigma }}{\tilde{a}}^H\left({\hat{\theta }}_p^{\left(k\right)}\right)R^{-1}\tilde{a}\left({\hat{\theta }}_p^{\left(k\right)}\right).$

7.3判断是否收敛

$J_k={\left({\hat{\delta }}^{(k+1)}\right)}^H{Q}_k{\hat{\delta }}^{(k+1)}$

在满足J_k-1-J_k>ε(ε为预设的门限值)时,继续迭代；否则,结束迭代。最后一次迭代 得到的角度估计值和幅相误差估计值为

步骤8:在距离维上进行迭代，固定目标角度为即为 步骤7中估计所得到的角度估计值，则可以在距离维上进行一维幅相误差校正。进行与步骤 7相同的迭代过程，可求得距离估计值和幅相误差估计值

同样定义代价函数

$J_k=\underset{p=1}{\overset{P}{\Sigma }}{a}^H({\hat{\theta }}^{\left(k\right)},{{\hat{r}}_p}^{\left(k\right)}){\left({{\Gamma }_R}^{\left(k\right)}\right)}^H{R}^{-1}{{\Gamma }_r}^{\left(k\right)}a({\hat{\theta }}^{\left(k\right)},{{\hat{r}}_p}^{\left(k\right)})$

式中：“||·||”表示Frobenius范数。

8.1定义空间谱公式为:

$p\left({\hat{r}}^{\left(k\right)}|{{\hat{\Gamma }}_R}^{\left(k\right)},{\hat{\theta }}^{\left(k\right)}\right)=\frac{1}{{\hat{a}}^H\left({\hat{\theta }}^{\left(k\right)},{\hat{r}}^{\left(k\right)}\right){\left({{\hat{\Gamma }}_R}^{\left(k\right)}\right)}^H{R}^{-1}{{\hat{\Gamma }}_R}^{\left(k\right)}\hat{a}\left({\hat{\theta }}^{\left(k\right)},{\hat{r}}^{\left(k\right)}\right)}$

根据利用上述空间谱公式求取对应的P个最大值时的距离便可 得到第k次迭代的P个目标距离的估计值。

8.2令w＝[1,0,...,0]^T，则满足约束条件则利用拉格朗日乘子法， 可求得幅相误差矩阵如下:

$\left(\begin{array}{c}{{\hat{\Gamma }}_{\theta }}^{(k+1)}=diag\left\{{\hat{\delta }}^{(k+1)}\right\}\\ {\hat{\delta }}^{(k+1)}=\frac{{\left(Q_k\right)}^{-1}w}{w^T{\left(Q_k\right)}^{-1}w}\end{array}\right)$

其中，$Q_k=\underset{p=1}{\overset{P}{\Sigma }}{\tilde{a}}^H\left({\hat{R}}_p^{\left(k\right)}\right)R^{-1}\tilde{a}\left({\hat{R}}_p^{\left(k\right)}\right).$

8.3判断是否收敛

$J_k={\left({\hat{\delta }}^{(k+1)}\right)}^H{Q}_k{\hat{\delta }}^{(k+1)}$

在满足J_k-1-J_k>ε(ε为预设的门限值)时,继续迭代；否则,结束迭代。最后一次迭代 得到的角度估计值和幅相误差估计值为

步骤9

由于阵列幅相误差来自阵列通道的不一致性，比较两个维度上估计出的幅相误差 对角矩阵，若满足则停止迭代，即得到幅相误差的估计值否则，令k＝ k+1，转到步骤7，继续迭代。

本发明的效果可以通过下列仿真实验加以说明：

仿真条件

本例中的无源定位阵列采用步骤2中所述方法产生的稀布阵列，如图1所示，沿着x 轴方向，在0～100m范围内稀疏布置8个阵元，近场目标窄带不相关信号源个数P＝1，位置分 别(-22.5°,1700m)，波长λ＝0.3m，信噪比均为SNR＝20dB，快拍数L＝400，角度维在-45°～ 45°范围内以0.1°为步长进行搜索，距离维在1200m～2000m范围内以1m为步长进行搜索。设 相位均方根误差δ_p＝10°，幅度均方根误差δ_a＝1dB，即

$\left(\begin{array}{cc}\Delta \phi ~U(-\sqrt{12{\sigma }_{\phi }^2}/2,\sqrt{12{\sigma }_{\phi }^2}/2)& {\sigma }_{\phi }={\delta }_{\phi }\frac{\pi }{180}\end{array}\right)$

$\left(\begin{array}{cc}\Delta a~U(-\sqrt{12{\sigma }_a^2}/2,\sqrt{12{\sigma }_a^2}/2)& {\delta }_a=10{\mathrm{lg\sigma }}_a\end{array}\right)$

计算可得相位误差△φ～U(-0.3023，0.3023)，单位rad，或者△φ～U(-17.3°, 17.3°)幅度误差△a～U(-2.1085，2.1085)，即接收信号幅度为a，则加入幅度误差后的幅度 a'＝a+△a。

仿真结果

图3给出了存在幅相误差的情况下，近场信号无源定位的二维方向图，旁瓣峰值达 到-5dB左右。采用本发明的方法进行幅相误差校正之后，获得的二维方向图如4所示，旁瓣 峰值降低到-20dB左右，有非常明显的校正效果。图5和图6还分别给出了，校正误差前后在 二维方向图的峰值位置上距离维和角度维的切面图对比，显然，经过本发明的幅相误差校 正后，即使在存在幅相误差的情况下仍然有效保证了近场无源定位二维方向图的性能。