CN104199021A - 基于调频广播信号的外辐射源雷达修正比相测角方法 - Google Patents

Abstract

本发明属于雷达目标方位角测量技术领域，特别涉及基于调频广播信号的外辐射源雷达修正比相测角方法。该基于调频广播信号的外辐射源雷达修正比相测角方法包括以下步骤：利用四单元均匀圆阵Adcock天线来接收外辐射源调频广播信号，利用参考天线接收参考信号；利用参考信号构建杂波空间，根据杂波空间，对外辐射源调频广播信号和参考信号，进行杂波相消处理，得出每个单元均匀圆阵Adcock天线对应的杂波相消后信号；对第k单元均匀圆阵Adcock天线对应的杂波相消后信号进行距离多普勒处理，根据距离多普勒处理结果，建立关于目标方位角的优化模型，求解所述关于目标方位角的优化模型，得出目标方位角。

Description

基于调频广播信号的外辐射源雷达修正比相测角方法

技术领域

本发明属于雷达目标方位角测量技术领域，特别涉及基于调频广播信号的外辐射源雷达修正比相测角方法。

背景技术

对于外辐射源雷达来说，利用外辐射源发射电磁波，而目标本身不直接辐射电磁能量。在外辐射源雷达工作时，通过天线接收来自外部的非协同辐射源，即第三方的直射波，以及该外部辐射源照射目标后形成的反射波或散射波，利用其携带的多普勒频移、多站接收信号的时间差和到达角等信息，提取出目标信息并消除无用信息和干扰，从而完成对目标定位。

在现有的利用外辐射源雷达实现雷达目标方位角估计的技术方案中，当参考信号中含有目标回波时，线测角方法的性能变得很差。这是由于当参考信号中含有目标回波时，在杂波相消过程中，参考信号中的目标回波将混入相消结果中，与天线接收信号中的目标回波共同参与后续的距离多普勒处理过程，将极大程度的影响距离多普勒处理结果的相位。

发明内容

本发明的目的在于提出基于调频广播信号的外辐射源雷达修正比相测角方法。在参考信号中含有目标回波的情况下，采用本发明能够精确地测量得到目标的方位角，提高目标定位精度，满足实际应用的要求。

本发明的实现思路是，采用均匀四圆阵Adcock天线对目标进行测角是有很多优点，但是当参考信号中含有目标回波时，这种测角精度迅速恶化。所以本发明在基于采用均匀四圆阵Adcock天线的基础上，对测角方法进行改进，来满足测角精度。

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

基于调频广播信号的外辐射源雷达修正比相测角方法包括以下步骤：

步骤1，利用四单元均匀圆阵Adcock天线来接收外辐射源调频广播信号，所述四单元均匀圆阵Adcock天线包括第1单元均匀圆阵Adcock天线至第4单元均匀圆阵Adcock天线，第k单元均匀圆阵Adcock天线接收的外辐射源调频广播信号为S_k，k取1至4，S_k为N维列向量，N为第m单元均匀圆阵Adcock天线接收的信号的长度；则四单元均匀圆阵Adcock天线来接收外辐射源调频广播信号S为：S＝[S₁,S₂,S₃,S₄]^T，上标T表示矩阵或向量的转置；

利用参考天线接收参考信号，参考天线接收的参考信号为S_ref，S_ref为N维行向量；

步骤2，利用参考信号S_ref构建杂波空间r，杂波空间r为N×(P+1)维的矩阵，P为设定的杂波对消阶数；

步骤3，根据构建的杂波空间r，对每个单元均匀圆阵Adcock天线接收的外辐射源调频广播信号和参考天线接收的参考信号，进行杂波相消处理，得出每个单元均匀圆阵Adcock天线对应的杂波相消后信号e_k；

步骤4，对第k单元均匀圆阵Adcock天线对应的杂波相消后信号e_k进行距离多普勒处理，得到对应的距离多普勒矩阵R-D_k，k取1至4；得出矩阵R-D_k中每个元素对应的幅度值，矩阵R-D_k中各个元素对应的幅度值的最大值为矩阵R-D_k的距离多普勒平面峰值，矩阵R-D_k的距离多普勒平面峰值对应的元素位于矩阵R-D_k的第τ_k行第f_d,k列，将矩阵R-D_k第τ_k行第f_d,k列的元素表示为R-D_k(τ_k,f_d,k)；构造矩阵4×1维的矩阵R以及4×1维的向量b：

$R={\left[\begin{array}{c}R-D_1({\mathrm{\τ}}_1,f_{d,1})\\ R-D_2({\mathrm{\τ}}_2,f_{d,2})\\ R-D_3({\mathrm{\τ}}_3,f_{d,3})\\ R-D_4({\mathrm{\τ}}_4,f_{d,4})\end{array}\right]}_{4\×1}$

$b={\left[\begin{array}{c}{\mathrm{\β}}^1\left(1\right)\\ {\mathrm{\β}}^2\left(1\right)\\ {\mathrm{\β}}^3\left(1\right)\\ {\mathrm{\β}}^4\left(1\right)\end{array}\right]}_{4\×1}$

其中，β^k(1)表示向量β^k的第1个元素，β^k＝(r^Tr)^-1r^TS_k，上标T表示矩阵或向量的转置，上标-1表示矩阵的逆；

步骤5，根据步骤4得出的矩阵R和向量b，建立关于目标方位角的优化模型，求解所述关于目标方位角的优化模型，得出目标方位角。

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

在所述步骤2中，参考信号S_ref为：

S_ref＝[S_ref(1),S_ref(2),...,S_ref(N)]

其中，S_ref(n)表示参考信号S_ref的第n列，n取1至N；

利用参考信号S_ref构建的杂波空间r为：

$r={\left[\begin{array}{cccc}& & & 0\\ & & & \·\\ & 0& & \·\\ & & & \·\\ S_{\mathrm{ref}}{(1~N)}^T& S_{\mathrm{ref}}{(1~N-1)}^T& \·\·\·& 0\\ & & & S_{\mathrm{ref}}{(1~N-P)}^T\end{array}\right]}_{N\×(P+1)}$

其中，S_ref(1～p)＝[S_ref(1),...,S_ref(p)]，p取N至N-P；上标T表示矩阵或向量的转置；所述杂波空间r的第1列为S_ref(1～N)^T，设p'为自然数且p'取2至P+1，杂波空间r的第p'列的前p'-1个元素均为零，杂波空间r的第p'列的后N-p'+1个元素依次为S_ref(1)至S_ref(N-p'+1)。

在步骤3中，第k单元均匀圆阵Adcock天线对应的杂波相消后信号e_k为：

e_k＝S_k-rβ^k

其中，k取1至4，β^k＝(r^Tr)^-1r^TS_k，上标T表示矩阵或向量的转置，上标-1表示矩阵的逆。

所述步骤4的具体子步骤为：

(4.1)确定最大检测距离单元个数NumRange和最大检测多普勒单元个数NumDoppler，NumRange<N；

(4.2)构造每个单元均匀圆阵Adcock天线对应的距离多普勒矩阵，第k单元均匀圆阵Adcock天线对应的距离多普勒矩阵为：

其中，k取1至4，R-D_k表示第k单元均匀圆阵Adcock天线对应的距离多普勒矩阵，矩阵R-D_k为NumDoppler×NumRange维的全零矩阵；

(4.3)设置变量i，i＝0,1,2,...,NumRange-1，当i＝0时，执行子步骤(4.4)；

(4.4)构造N维行向量stemp，所述N维行向量stemp的前i个元素均为零，N维行向量stemp的后N-i个元素依次为S_ref(1)至S_ref(N-i)；S_ref(n)表示参考信号S_ref的第n列，n取1至N；

(4.5)构造NumDoppler维的行向量

${\mathrm{res}}_1^k=[{\mathrm{res}}^k\left(1\right),{\mathrm{res}}^k(1+M),{\mathrm{res}}^k(1+2M),\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{res}}^k(1+(\mathrm{NumDoppler}-1)\&times;M)]$

其中，k取1至4，res^k(j)＝stemp(j)×e_k(j)，j＝1,2,…,N，stemp(j)表示行向量stemp的第j个元素，e_k(j)表示第k单元均匀圆阵Adcock天线对应的杂波相消后信号e_k的第j个元素；M表示设置的抽头系数，M满足1+(NumDoppler-1)×M≤N；

(4.6)对行向量作FFT变换，得到向量对进行倒置得到列向量然后更新第k单元均匀圆阵Adcock天线对应的距离多普勒矩阵R-D_k；所述更新矩阵R-D_k的过程为：将矩阵R-D_k的第i+1列用列向量代替；

(4.7)判断i是否等于NumRange，如果是，则得出经更新后的第k单元均匀圆阵Adcock天线对应的距离多普勒矩阵R-D_k，执行子步骤(4.8)，如果否，则令i的值自增1，返回至子步骤(4.4)；

(4.8)将k从1至4进行遍历，重复执行子步骤(4.3)至子步骤(4.7)，得出经更新后的每个单元均匀圆阵Adcock天线对应的R-D矩阵；

(4.9)得出矩阵R-D_k中每个元素对应的幅度值，矩阵R-D_k中各个元素对应的幅度值的最大值为矩阵R-D_k的距离多普勒平面峰值，矩阵R-D_k的距离多普勒平面峰值对应的元素位于矩阵R-D_k的第τ_k行第f_d,k列，将矩阵R-D_k第τ_k行第f_d,k列的元素表示为R-D_k(τ_k,f_d,k)；

构造矩阵4×1维的矩阵R以及4×1维的向量b：

$R={\left[\begin{array}{c}R-D_1({\mathrm{\&tau;}}_1,f_{d,1})\\ R-D_2({\mathrm{\&tau;}}_2,f_{d,2})\\ R-D_3({\mathrm{\&tau;}}_3,f_{d,3})\\ R-D_4({\mathrm{\&tau;}}_4,f_{d,4})\end{array}\right]}_{4\&times;1}$

$b={\left[\begin{array}{c}{\mathrm{\&beta;}}^1\left(1\right)\\ {\mathrm{\&beta;}}^2\left(1\right)\\ {\mathrm{\&beta;}}^3\left(1\right)\\ {\mathrm{\&beta;}}^4\left(1\right)\end{array}\right]}_{4\&times;1}$

其中，β^k(1)表示向量β^k的第1个元素。

所述步骤5的具体子步骤为：

(5.1)建立如下关于目标方位角的优化模型：

${\mathrm{\&theta;}}_s=\underset{\mathrm{\&theta;}}{\min }\left|\right|R+\mathrm{b\&xi;}-\mathrm{Ta}\left(\mathrm{\&theta;}\right)\left|\right|=\underset{\mathrm{\&theta;}}{\min }\mathrm{\&psi;}(\mathrm{\&theta;},T,\mathrm{\&xi;})$

其中，ψ(θ,T,ξ)＝||R+bξ-Ta(θ)||，||·||表示向量的模，ξ、T均为数值变量，θ为角度变量，θ_s表示目标的方位角；a(θ)为：

$a\left(\mathrm{\&theta;}\right)=\left[\begin{array}{c}{e}^{j\frac{2\mathrm{\&pi;l}}{\mathrm{\&lambda;}}\cos \mathrm{\&theta;}}\\ e^{j\frac{2\mathrm{\&pi;l}}{\mathrm{\&lambda;}}\cos (\mathrm{\&theta;}-\frac{\mathrm{\&pi;}}{2})}\\ e^{j\frac{2\mathrm{\&pi;l}}{\mathrm{\&lambda;}}\cos (\mathrm{\&theta;}-\mathrm{\&pi;})}\\ e^{j\frac{2\mathrm{\&pi;l}}{\mathrm{\&lambda;}}\cos (\mathrm{\&theta;}-\frac{3\mathrm{\&pi;}}{2})}\end{array}\right]$

其中，l为每个单元均匀圆阵Adcock天线的圆阵半径，λ为外辐射源的调频广播信号的波长；

(5.2)对子步骤(5.1)中建立的目标方位角的优化模型进行求解，得出目标方位角θ_s。

在子步骤(5.2)中，对子步骤(5.1)中建立的目标方位角的优化模型进行求解的过程包括以下子步骤：

a)根据公式 $\left[\begin{array}{c}{T}^{*}\\ {\mathrm{\&xi;}}^{*}\end{array}\right]={\left(A^H\left(\mathrm{\&theta;}\right)A\left(\mathrm{\&theta;}\right)\right)}^{-1}{A}^H\left(\mathrm{\&theta;}\right)R,$ 计算数值变量T与ξ的最小二乘解，其中，上标*表示共轭，上标H表示矩阵的共轭转置，上标-1表示矩阵的逆，A(θ)＝[a(θ)-b]；将求解出的T与ξ的最小二乘解代入ψ(θ,T,ξ)中，得出函数ψ(θ)；

b)设定角度初始值θ₁，设定允许误差

c)设置参数H₁＝1，计算ψ(θ)在角度初始值θ₁处的梯度g₁，设置整数变量m，m＝1,2...，当m＝1时，执行子步骤d)；

d)计算搜索方向d_m，d_m＝-H_mg_m；

e)求解迭代步长λ_m，使其满足 $\mathrm{\&psi;}({\mathrm{\&theta;}}_m+{\mathrm{\&lambda;}}_m{d}_m)=\underset{{\mathrm{\&lambda;}}^{\&prime;}\&GreaterEqual;0}{\min }\mathrm{\&psi;}({\mathrm{\&theta;}}_m+{\mathrm{\&lambda;}}^{\&prime;}{d}_m),$ 其中，λ'为大于等于0的变量，表示：当λ'在[0,+∞)变化时ψ(θ_m+λ'd_m)能取到的最大值；然后，令θ_m+1＝θ_m+λ_md_m；

f)判断与允许误差的大小关系，表示ψ(θ)在θ_m+1处的梯度，表示的绝对值；若则得出θ_s＝θ_m+1；否则进行子步骤g)；

g)令p_m＝θ_m+1-θ_m，q_m＝g_m+1-g_m，然后计算H_m+1，

$H_{m+1}=H_m+\frac{p_m{p}_m}{p_m{q}_m}-\frac{H_m{q}_m{q}_m{H}_m}{q_m{H}_m{q}_m}$

然后，令m的值自增1，返回子步骤d)。

本发明的有益效果为：1)本发明在采用均匀四圆阵Adcock天线的基础上，对测角方法进行改进，所以该发明有Adcock天线架构小，便于安装和移动，天线的阵子个数少，阵子之间的互耦和效应相对较弱和测量角度不模糊的优点。2)参考信号中含有目标回波，当参考信号的直信比发生变化时，天线接收信号中目标回波的信噪比发生变化时，目标的方位角发生变化时，本发明均能较精确地测量目标方位角，而参考信号中目标回波的强度大到一定程度时，原始方法的估计性能急剧恶化，其应用范围大大受限，所以本发明应用范围更加广泛。3)当零延迟附近存在强多径信号时，采用本发明的修正测角方法，使得目标测角几乎不受强多径信号多径延迟大小的影响。

附图说明

图1为本发明的基于调频广播信号的外辐射源雷达修正比相测角方法的流程图；

图2为仿真1中两种方法得出的方位角估计均方根误差随着参考信号中直信比的变化关系示意图；

图3a是仿真2中本发明的方位角估计均方根误差随天线接收信号中目标信噪比的变化关系示意图；

图3b为仿真2中传统方法方位角估计均方根误差随天线接收信号中目标回波信噪比的变化关系示意图。

具体实施方式

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

参照图1，为本发明的基于调频广播信号的外辐射源雷达修正比相测角方法的流程图。该基于调频广播信号的外辐射源雷达修正比相测角方法包括以下步骤：

步骤1，利用四单元均匀圆阵Adcock天线来接收外辐射源调频广播信号，所述外辐射源调频广播信号包括：外辐射源的直达波信号(外辐射源直接发送至四单元均匀圆阵Adcock天线的调频广播信号)、外辐射源的多径信号(外辐射源通过多径发送至四单元均匀圆阵Adcock天线的调频广播信号)、以及外辐射源的目标回波信号(外辐射源的调频广播信号照射到目标后被四单元均匀圆阵Adcock天线接收的目标回波信号)。四单元均匀圆阵Adcock天线包括第1单元均匀圆阵Adcock天线至第4单元均匀圆阵Adcock天线，第k单元均匀圆阵Adcock天线接收的外辐射源调频广播信号为S_k，k取1至4，S_k为N维列向量，N为第m单元均匀圆阵Adcock天线接收的信号的长度；则四单元均匀圆阵Adcock天线来接收外辐射源调频广播信号S为：S＝[S₁,S₂,S₃,S₄]^T，上标T表示矩阵或向量的转置。

利用参考天线接收含目标回波的参考信号，参考天线接收的参考信号为S_ref，S_ref为N维行向量。

步骤2，利用参考信号S_ref构建杂波空间r，杂波空间r为N×(P+1)维的矩阵，P为设定的杂波对消阶数。

具体地，参考信号S_ref为：

S_ref＝[S_ref(1),S_ref(2),...,S_ref(N)]

其中，S_ref(n)表示参考信号S_ref的第n列，n取1至N。

利用参考信号S_ref构建的杂波空间r为：

$r={\left[\begin{array}{cccc}& & & 0\\ & & & \&CenterDot;\\ & 0& & \&CenterDot;\\ & & & \&CenterDot;\\ S_{\mathrm{ref}}{(1~N)}^T& S_{\mathrm{ref}}{(1~N-1)}^T& \&CenterDot;\&CenterDot;\&CenterDot;& 0\\ & & & S_{\mathrm{ref}}{(1~N-P)}^T\end{array}\right]}_{N\&times;(P+1)}$

其中，S_ref(1～p)表示由参考信号S_ref的第1列至第p列组成的向量，p取N至N-P，即S_ref(1～p)＝[S_ref(1),...,S_ref(p)]；上标T表示矩阵或向量的转置；所述杂波空间r的第1列为S_ref(1～N)^T，设p'为自然数且p'取2至P+1，则杂波空间r的第p'列的前p'-1个元素均为零，杂波空间r的第p'列的后N-p'+1个元素依次为S_ref(1)至S_ref(N-p'+1)。

本发明实施例中，将P经验性地设置为256、144或64。

步骤3，根据构建的杂波空间r，对每个单元均匀圆阵Adcock天线接收的外辐射源调频广播信号和参考天线接收的参考信号，进行杂波相消处理，得出每个单元均匀圆阵Adcock天线对应的杂波相消后信号。

具体地，第k单元均匀圆阵Adcock天线对应的杂波相消后信号e_k为：

e_k＝S_k-rβ^k

其中，k取1至4，β^k＝(r^Tr)^-1r^TS_k，上标T表示矩阵或向量的转置，上标-1表示矩阵的逆。

步骤4，对第k单元均匀圆阵Adcock天线对应的杂波相消后信号e_k进行距离多普勒处理，得到对应的距离多普勒矩阵R-D_k，k取1至4。得出矩阵R-D_k中每个元素对应的幅度值，矩阵R-D_k中各个元素对应的幅度值的最大值为矩阵R-D_k的距离多普勒平面峰值，矩阵R-D_k的距离多普勒平面峰值对应的元素位于矩阵R-D_k的第τ_k行第f_d,k列，将矩阵R-D_k第τ_k行第f_d,k列的元素表示为R-D_k(τ_k,f_d,k)。构造矩阵4×1维的矩阵R以及4×1维的向量b：

$R={\left[\begin{array}{c}R-D_1({\mathrm{\&tau;}}_1,f_{d,1})\\ R-D_2({\mathrm{\&tau;}}_2,f_{d,2})\\ R-D_3({\mathrm{\&tau;}}_3,f_{d,3})\\ R-D_4({\mathrm{\&tau;}}_4,f_{d,4})\end{array}\right]}_{4\&times;1}$

$b={\left[\begin{array}{c}{\mathrm{\&beta;}}^1\left(1\right)\\ {\mathrm{\&beta;}}^2\left(1\right)\\ {\mathrm{\&beta;}}^3\left(1\right)\\ {\mathrm{\&beta;}}^4\left(1\right)\end{array}\right]}_{4\&times;1}$

其中，β^k(1)表示向量β^k的第1个元素，β^k＝(r^Tr)^-1r^TS_k，上标T表示矩阵或向量的转置，上标-1表示矩阵的逆。

其具体子步骤为：

(4.1)确定最大检测距离单元个数NumRange和最大检测多普勒单元个数NumDoppler，NumRange<N。本发明实施例中，将最大检测距离单元个数NumRange经验性地设置为256，将最大检测多普勒单元个数NumDoppler经验性地设置为2048。

(4.2)构造每个单元均匀圆阵Adcock天线对应的距离多普勒矩阵，第k单元均匀圆阵Adcock天线对应的距离多普勒矩阵为：

其中，k取1至4，R-D_k表示第k单元均匀圆阵Adcock天线对应的距离多普勒矩阵，矩阵R-D_k为NumDoppler×NumRange维的全零矩阵；

(4.3)设置变量i，i＝0,1,2,...,NumRange-1，当i＝0时，执行子步骤(4.4)；

(4.4)构造N维行向量stemp：

$\mathrm{stemp}={\left[\begin{array}{ccccccc}0& 0& ...& 0& S_{\mathrm{ref}}\left(1\right)& S_{\mathrm{ref}}\left(2\right)& S_{\mathrm{ref}}(N-i)\end{array}\right]}_{1\&times;N}$

N维行向量stemp的前i个元素均为零，N维行向量stemp的后N-i个元素依次为S_ref(1)至S_ref(N-i)。

(4.5)构造NumDoppler维的行向量

${\mathrm{res}}_1^k=[{\mathrm{res}}^k\left(1\right),{\mathrm{res}}^k(1+M),{\mathrm{res}}^k(1+2M),\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{res}}^k(1+(\mathrm{NumDoppler}-1)\&times;M)]$

其中，k取1至4，res^k(j)＝stemp(j)×e_k(j)，j＝1,2,…,N，stemp(j)表示行向量stemp的第j个元素，e_k(j)表示第k单元均匀圆阵Adcock天线对应的杂波相消后信号e_k的第j个元素；M表示设置的抽头系数，M满足1+(NumDoppler-1)×M≤N的条件，本发明实施例中，经验性地将M设置为100。

(4.6)对行向量作FFT变换，得到向量 FFT(·)表示进行FFT变换。对进行倒置得到列向量然后更新第k单元均匀圆阵Adcock天线对应的距离多普勒矩阵R-D_k；所述更新矩阵R-D_k的过程为：将矩阵R-D_k的第i+1列用列向量代替。

(4.7)判断i是否等于NumRange，如果是，则得出经更新后的第k单元均匀圆阵Adcock天线对应的距离多普勒矩阵R-D_k，执行子步骤(4.8)，如果否，则令i的值自增1，返回至子步骤(4.4)。

(4.8)将k从1至4进行遍历，重复执行子步骤(4.3)至子步骤(4.7)，得出经更新后的每个单元均匀圆阵Adcock天线对应的R-D矩阵。

(4.9)对矩阵R-D_k进行恒虚警检测，得出矩阵R-D_k的距离多普勒平面峰值，根据矩阵R-D_k的距离多普勒平面峰值对应的距离多普勒单元(τ_k,f_d,k)，得出矩阵R-D_k第τ_k行第f_d,k列的元素R-D_k(τ_k,f_d,k)，k取1至4。

具体地说，首先得出矩阵R-D_k中每个元素对应的幅度值(矩阵R-D_k中每个元素的值)，矩阵R-D_k中各个元素对应的幅度值的最大值为矩阵R-D_k的距离多普勒平面峰值，矩阵R-D_k的距离多普勒平面峰值对应的元素位于矩阵R-D_k的第τ_k行第f_d,k列，将矩阵R-D_k第τ_k行第f_d,k列的元素表示为R-D_k(τ_k,f_d,k)。

然后，构造矩阵4×1维的矩阵R以及4×1维的向量b：

$R={\left[\begin{array}{c}R-D_1({\mathrm{\&tau;}}_1,f_{d,1})\\ R-D_2({\mathrm{\&tau;}}_2,f_{d,2})\\ R-D_3({\mathrm{\&tau;}}_3,f_{d,3})\\ R-D_4({\mathrm{\&tau;}}_4,f_{d,4})\end{array}\right]}_{4\&times;1}$

$b={\left[\begin{array}{c}{\mathrm{\&beta;}}^1\left(1\right)\\ {\mathrm{\&beta;}}^2\left(1\right)\\ {\mathrm{\&beta;}}^3\left(1\right)\\ {\mathrm{\&beta;}}^4\left(1\right)\end{array}\right]}_{4\&times;1}$

其中，β^k(1)表示向量β^k的第1个元素。

步骤5，根据距离多普勒处理结果，建立关于目标方位角的优化模型，求解所述关于目标方位角的优化模型，得出目标方位角。

其具体子步骤为：

(5.1)建立如下关于目标方位角的优化模型：

${\mathrm{\&theta;}}_s=\underset{\mathrm{\&theta;}}{\min }\left|\right|R+\mathrm{b\&xi;}-\mathrm{Ta}\left(\mathrm{\&theta;}\right)\left|\right|=\underset{\mathrm{\&theta;}}{\min }\mathrm{\&psi;}(\mathrm{\&theta;},T,\mathrm{\&xi;})$

其中，ψ(θ,T,ξ)＝||R+bξ-Ta(θ)||，||·||表示向量的模，ξ、T均为数值变量，θ为角度变量，θ_s表示目标的方位角。a(θ)为导向矢量，a(θ)为：

$a\left(\mathrm{\&theta;}\right)=\left[\begin{array}{c}{e}^{j\frac{2\mathrm{\&pi;l}}{\mathrm{\&lambda;}}\cos \mathrm{\&theta;}}\\ e^{j\frac{2\mathrm{\&pi;l}}{\mathrm{\&lambda;}}\cos (\mathrm{\&theta;}-\frac{\mathrm{\&pi;}}{2})}\\ e^{j\frac{2\mathrm{\&pi;l}}{\mathrm{\&lambda;}}\cos (\mathrm{\&theta;}-\mathrm{\&pi;})}\\ e^{j\frac{2\mathrm{\&pi;l}}{\mathrm{\&lambda;}}\cos (\mathrm{\&theta;}-\frac{3\mathrm{\&pi;}}{2})}\end{array}\right]$

其中，l为每个单元均匀圆阵Adcock天线的圆阵半径，λ为外辐射源的调频广播信号的波长。

上述优化模型的求解过程可以看作：通过求解非线性方程

$\left\{\begin{array}{c}R\left(1\right)={\mathrm{Ta}}^1\left(\mathrm{\&theta;}\right)-b\left(1\right)\mathrm{\&xi;}\\ R\left(2\right)={\mathrm{Ta}}^2\left(\mathrm{\&theta;}\right)-b\left(2\right)\mathrm{\&xi;}\\ R\left(3\right)={\mathrm{Ta}}^3\left(\mathrm{\&theta;}\right)-b\left(3\right)\mathrm{\&xi;}\\ R\left(4\right)={\mathrm{Ta}}^4\left(\mathrm{\&theta;}\right)-b\left(4\right)\mathrm{\&xi;}\end{array}\right.$

来得出目标的方位角θ_s。其中，R(k)表示矩阵R的第k行，k取1至4，b(k)表示向量b的第k行，a^k(θ)表示向量a(θ)的第k个元素。

(5.2)对子步骤(5.1)中建立的目标方位角的优化模型进行求解，得出目标方位角θ_s。

对子步骤(5.1)中建立的目标方位角的优化模型进行求解的过程包括以下子步骤：

a)根据公式 $\left[\begin{array}{c}{T}^{*}\\ {\mathrm{\&xi;}}^{*}\end{array}\right]={\left(A^H\left(\mathrm{\&theta;}\right)A\left(\mathrm{\&theta;}\right)\right)}^{-1}{A}^H\left(\mathrm{\&theta;}\right)R,$ 计算数值变量T与ξ的最小二乘解，其中，上标*表示共轭，上标H表示矩阵的共轭转置，上标-1表示矩阵的逆，A(θ)＝[a(θ)-b]。将求解出的T与ξ的最小二乘解代入ψ(θ,T,ξ)中，得出函数ψ(θ)。

b)设定角度初始值θ₁，设定允许误差

c)设置参数H₁＝1，计算ψ(θ)在角度初始值θ₁处的梯度g₁，设置整数变量m，m＝1,2...，当m＝1时，执行子步骤d)。

d)计算搜索方向d_m，d_m＝-H_mg_m。

e)从θ_m出发，沿着d_m方向搜索，求迭代步长λ_m，使其满足 $\mathrm{\&psi;}({\mathrm{\&theta;}}_m+{\mathrm{\&lambda;}}_m{d}_m)=\underset{{\mathrm{\&lambda;}}^{\&prime;}\&GreaterEqual;0}{\min }\mathrm{\&psi;}({\mathrm{\&theta;}}_m+{\mathrm{\&lambda;}}^{\&prime;}{d}_m),$ 其中，λ'为大于等于0的变量，表示：当λ'在[0,+∞)变化时ψ(θ_m+λ'd_m)能取到的最大值；然后，令θ_m+1＝θ_m+λ_md_m；

f)判断与允许误差的大小关系，表示ψ(θ)在θ_m+1处的梯度，表示的绝对值；若则停止迭代，得出θ_s＝θ_m+1；否则进行子步骤g)；

g)令p_m＝θ_m+1-θ_m，q_m＝g_m+1-g_m，然后计算H_m+1，

$H_{m+1}=H_m+\frac{p_m{p}_m}{p_m{q}_m}-\frac{H_m{q}_m{q}_m{H}_m}{q_m{H}_m{q}_m}$

然后，令m的值自增1，返回子步骤d)。

本发明的效果可通过以下仿真实验进一步说明：

1)仿真实验条件：

在实验中，硬件平台选用HP Z820工作站，Win7系统，软件平台为MATLAB2010。仿真实验采用一帧外辐射源雷达仿真数据，其中主通道中包括目标信号、直达波信号、多径信号和噪声信号。天线，直达波，多径杂波和目标回波的信息设置如实验效果所示。在仿真实验中，采用本发明和传统方法对目标方位角进行估计。传统方法指在专利《基于Adcock天线的外辐射源相干定位雷达测向系统及其测量方法》(公开号CN101799542A，专利号201010013570，申请日2010.01.08)中的目标方位角估计方法。

2)实验效果：

仿真1：方位角估计均方根误差随参考信号直信比的变化对比图。在仿真1中，设定天线接收信号中直达波的信噪比为40dB，直达波的方位角为0°；多径杂波的信噪比为10dB，多径延迟为3个距离单元，方位角为30°；目标回波的信噪比为-10dB，直信比为50dB，方位角为60°，目标回波所在距离多普勒单元为(60,102)。参考信号中直达波与目标回波的能量比(直信比)从50dB逐渐增大到80dB，即参考信号的直信比比天线接收信号的直信比大0到30dB。目标方位角均方根误差的计算公式为：

${\widehat{\mathrm{\&theta;}}}_s=\sqrt{\frac{1}{N^{\&prime;}}\underset{i=1}{\overset{N^{\&prime;}}{\mathrm{\&Sigma;}}}{({\widetilde{\mathrm{\&theta;}}}_s^i-{\mathrm{\&theta;}}_s)}^2}$

其中，N'为蒙特卡洛仿真次数，本文仿真实验中该值均为100，为第i次蒙特卡洛仿真测得的目标方位角，θ_s为目标方位角真值。

根据上述条件得到的方位角估计均方根误差随着参考信号中直信比的变化关系如图2所示，图2中，横轴表示参考信号的直信比，单位为dB，纵轴表示目标方位角均方根误差，单位为度。以十字标示的曲线为传统方法中目标方位角均方根误差随参考信号直信比的变化曲线，以圆圈标示的曲线为本发明中目标方位角均方根误差随参考信号直信比的变化曲线。由图2可知，本发明的目标方位角均方根误差几乎不随参考信号的直信比变化而变化。当参考信号的直信比在50dB到70dB之间，即比天线接收信号的直信比大0到20dB时，传统方法的方位角估计精度极大程度地恶化。天线接收信号的直信比确定后，参考信号的直信比越高，传统方法的方位角估计精度就越高，随着直信比的下降，传统方法方位角估计精度迅速恶化。

仿真2：对不同目标回波信噪比时原始方法与改进方法的方位角测定进行对比，将参考信号的直信比设置为60dB，天线接收信号中目标回波的信噪比从-20dB逐渐增大到0dB，其他条件如仿真1，两种方法的方位角估计精度的变化曲线如图3所示。图3a是本发明的方位角估计均方根误差随天线接收信号中目标信噪比的变化关系示意图。图3a中，横轴表示天线接收信号中目标信噪比，单位为dB，纵轴表示目标方位角估计均方根误差，单位为度。如图3a所示，信噪比高时本发明方位角估计精度提高，信噪比降低时本发明的方位角估计精度降低。图3b为传统方法方位角估计均方根误差随天线接收信号中目标回波信噪比的变化关系示意图。图3b中，横轴表示天线接收信号中目标信噪比，单位为dB，纵轴表示目标方位角估计均方根误差，单位为度。由图3a和图3b的对比可知，对于每一个相同的目标信噪比，传统方法的方位角估计精度均远远差于本发明，天线接收信号中目标回波信噪比提高时，传统方法的均方根误差变小，趋近于本发明的精度。天线接收信号中目标回波的信噪比提高，传统方法方位角估计精度提高，这与噪声和对消的影响减弱也有关系。

仿真3：不同方位角时目标方位角估计均方根误差，将参考信号的直信比设置为60dB，天线接收信号中目标回波的信噪比设置为-10dB，其他条件如仿真1，此时参考信号的直信比比天线接收信号的直信比大10dB，不同目标方位角时的本发明和传统方法的方位角估计精度如表1所示，表1中所有数值的单位均为度，我们对表1中的每一个角度都进行了100次蒙特卡洛仿真，并计算对应的均方根误差。

表1

由表1可知，当目标的方位角发生变化时，改进方法总能精确测量得到目标的方位角。参考信号的直信比比天线接收信号的直信比大10dB时，原始方法的方位角估计性能远远差于改进方法。

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

Claims (6)

1.基于调频广播信号的外辐射源雷达修正比相测角方法，其特征在于，包括以下步骤：

步骤1，利用四单元均匀圆阵Adcock天线来接收外辐射源调频广播信号，所述四单元均匀圆阵Adcock天线包括第1单元均匀圆阵Adcock天线至第4单元均匀圆阵Adcock天线，第k单元均匀圆阵Adcock天线接收的外辐射源调频广播信号为S_k，k取1至4，S_k为N维列向量，N为第m单元均匀圆阵Adcock天线接收的信号的长度；则四单元均匀圆阵Adcock天线来接收外辐射源调频广播信号S为：S＝[S₁,S₂,S₃,S₄]^T，上标T表示矩阵或向量的转置；

利用参考天线接收参考信号，参考天线接收的参考信号为S_ref，S_ref为N维行向量；

步骤2，利用参考信号S_ref构建杂波空间r，杂波空间r为N×(P+1)维的矩阵，P为设定的杂波对消阶数；

步骤3，根据构建的杂波空间r，对每个单元均匀圆阵Adcock天线接收的外辐射源调频广播信号和参考天线接收的参考信号，进行杂波相消处理，得出每个单元均匀圆阵Adcock天线对应的杂波相消后信号；

步骤4，对第k单元均匀圆阵Adcock天线对应的杂波相消后信号e_k进行距离多普勒处理，得到对应的距离多普勒矩阵R-D_k，k取1至4；得出矩阵R-D_k中每个元素对应的幅度值，矩阵R-D_k中各个元素对应的幅度值的最大值为矩阵R-D_k的距离多普勒平面峰值，矩阵R-D_k的距离多普勒平面峰值对应的元素位于矩阵R-D_k的第τ_k行第f_d,k列，将矩阵R-D_k第τ_k行第f_d,k列的元素表示为R-D_k(τ_k,f_d,k)；构造矩阵4×1维的矩阵R以及4×1维的向量b：

$R={\left[\begin{array}{c}R-D_1({\mathrm{\&tau;}}_1,f_{d,1})\\ R-D_2({\mathrm{\&tau;}}_2,f_{d,2})\\ R-D_3({\mathrm{\&tau;}}_3,f_{d,3})\\ R-D_4({\mathrm{\&tau;}}_4,f_{d,4})\end{array}\right]}_{4\&times;1}$

$b={\left[\begin{array}{c}{\mathrm{\&beta;}}^1\left(1\right)\\ {\mathrm{\&beta;}}^2\left(1\right)\\ {\mathrm{\&beta;}}^3\left(1\right)\\ {\mathrm{\&beta;}}^4\left(1\right)\end{array}\right]}_{4\&times;1}$

其中，β^k(1)表示向量β^k的第1个元素，β^k＝(r^Tr)^-1r^TS_k，上标T表示矩阵或向量的转置，上标-1表示矩阵的逆；k取1至4；

步骤5，根据步骤4得出的矩阵R和向量b，建立关于目标方位角的优化模型，求解所述关于目标方位角的优化模型，得出目标方位角。

2.如权利要求1所述的基于调频广播信号的外辐射源雷达修正比相测角方法，其特征在于，在所述步骤2中，参考信号S_ref为：

S_ref＝[S_ref(1),S_ref(2),...,S_ref(N)]

其中，S_ref(n)表示参考信号S_ref的第n列，n取1至N；

利用参考信号S_ref构建的杂波空间r为：

$r={\left[\begin{array}{cccc}& & & 0\\ & & & \&CenterDot;\\ & 0& & \&CenterDot;\\ & & & \&CenterDot;\\ S_{\mathrm{ref}}{(1~N)}^T& S_{\mathrm{ref}}{(1~N-1)}^T& \&CenterDot;\&CenterDot;\&CenterDot;& 0\\ & & & S_{\mathrm{ref}}{(1~N-P)}^T\end{array}\right]}_{N\&times;(P+1)}$

其中，S_ref(1～p)＝[S_ref(1),...,S_ref(p)]，p取N至N-P；上标T表示矩阵或向量的转置；所述杂波空间r的第1列为S_ref(1～N)^T，设p'为自然数且p'取2至P+1，杂波空间r的第p'列的前p'-1个元素均为零，杂波空间r的第p'列的后N-p'+1个元素依次为S_ref(1)至S_ref(N-p'+1)。

3.如权利要求1所述的基于调频广播信号的外辐射源雷达修正比相测角方法，其特征在于，在步骤3中，第k单元均匀圆阵Adcock天线对应的杂波相消后信号e_k为：

e_k＝S_k-rβ^k

其中，k取1至4，β^k＝(r^Tr)^-1r^TS_k，上标T表示矩阵或向量的转置，上标-1表示矩阵的逆。

4.如权利要求1所述的基于调频广播信号的外辐射源雷达修正比相测角方法，其特征在于，所述步骤4的具体子步骤为：

(4.1)确定最大检测距离单元个数NumRange和最大检测多普勒单元个数NumDoppler，NumRange<N；

(4.2)构造每个单元均匀圆阵Adcock天线对应的距离多普勒矩阵，第k单元均匀圆阵Adcock天线对应的距离多普勒矩阵为：

其中，k取1至4，R-D_k表示第k单元均匀圆阵Adcock天线对应的距离多普勒矩阵，矩阵R-D_k为NumDoppler×NumRange维的全零矩阵；

(4.3)设置变量i，i＝0,1,2,...,NumRange-1，当i＝0时，执行子步骤(4.4)；

(4.4)构造N维行向量stemp，所述N维行向量stemp的前i个元素均为零，N维行向量stemp的后N-i个元素依次为S_ref(1)至S_ref(N-i)；S_ref(n)表示参考信号S_ref的第n列，n取1至N；

(4.5)构造NumDoppler维的行向量

${\mathrm{res}}_1^k=[{\mathrm{res}}^k\left(1\right),{\mathrm{res}}^k(1+M),{\mathrm{res}}^k(1+2M),\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{res}}^k(1+(\mathrm{NumDoppler}-1)\&times;M)]$

其中，k取1至4，res^k(j)＝stemp(j)×e_k(j)，j＝1,2,…,N，stemp(j)表示行向量stemp的第j个元素，e_k(j)表示第k单元均匀圆阵Adcock天线对应的杂波相消后信号e_k的第j个元素；M表示设置的抽头系数，M满足1+(NumDoppler-1)×M≤N；

(4.6)对行向量作FFT变换，得到向量对进行倒置得到列向量然后更新第k单元均匀圆阵Adcock天线对应的距离多普勒矩阵R-D_k；所述更新矩阵R-D_k的过程为：将矩阵R-D_k的第i+1列用列向量代替；

(4.7)判断i是否等于NumRange，如果是，则得出经更新后的第k单元均匀圆阵Adcock天线对应的距离多普勒矩阵R-D_k，执行子步骤(4.8)，如果否，则令i的值自增1，返回至子步骤(4.4)；

(4.8)将k从1至4进行遍历，重复执行子步骤(4.3)至子步骤(4.7)，得出经更新后的每个单元均匀圆阵Adcock天线对应的R-D矩阵；

(4.9)得出矩阵R-D_k中每个元素对应的幅度值，矩阵R-D_k中各个元素对应的幅度值的最大值为矩阵R-D_k的距离多普勒平面峰值，矩阵R-D_k的距离多普勒平面峰值对应的元素位于矩阵R-D_k的第τ_k行第f_d,k列，将矩阵R-D_k第τ_k行第f_d,k列的元素表示为R-D_k(τ_k,f_d,k)；

构造矩阵4×1维的矩阵R以及4×1维的向量b：

$R={\left[\begin{array}{c}R-D_1({\mathrm{\&tau;}}_1,f_{d,1})\\ R-D_2({\mathrm{\&tau;}}_2,f_{d,2})\\ R-D_3({\mathrm{\&tau;}}_3,f_{d,3})\\ R-D_4({\mathrm{\&tau;}}_4,f_{d,4})\end{array}\right]}_{4\&times;1}$

$b={\left[\begin{array}{c}{\mathrm{\&beta;}}^1\left(1\right)\\ {\mathrm{\&beta;}}^2\left(1\right)\\ {\mathrm{\&beta;}}^3\left(1\right)\\ {\mathrm{\&beta;}}^4\left(1\right)\end{array}\right]}_{4\&times;1}$

其中，β^k(1)表示向量β^k的第1个元素。

5.如权利要求1所述的基于调频广播信号的外辐射源雷达修正比相测角方法，其特征在于，所述步骤5的具体子步骤为：

(5.1)建立如下关于目标方位角的优化模型：

${\mathrm{\&theta;}}_s=\underset{\mathrm{\&theta;}}{\min }\left|\right|R+\mathrm{b\&xi;}-\mathrm{Ta}\left(\mathrm{\&theta;}\right)\left|\right|=\underset{\mathrm{\&theta;}}{\min }\mathrm{\&psi;}(\mathrm{\&theta;},T,\mathrm{\&xi;})$

其中，ψ(θ,T,ξ)＝||R+bξ-Ta(θ)||，||·||表示向量的模，ξ、T均为数值变量，θ为角度变量，θ_s表示目标的方位角；a(θ)为：

$a\left(\mathrm{\&theta;}\right)=\left[\begin{array}{c}{e}^{j\frac{2\mathrm{\&pi;l}}{\mathrm{\&lambda;}}\cos \mathrm{\&theta;}}\\ e^{j\frac{2\mathrm{\&pi;l}}{\mathrm{\&lambda;}}\cos (\mathrm{\&theta;}-\frac{\mathrm{\&pi;}}{2})}\\ e^{j\frac{2\mathrm{\&pi;l}}{\mathrm{\&lambda;}}\cos (\mathrm{\&theta;}-\mathrm{\&pi;})}\\ e^{j\frac{2\mathrm{\&pi;l}}{\mathrm{\&lambda;}}\cos (\mathrm{\&theta;}-\frac{3\mathrm{\&pi;}}{2})}\end{array}\right]$

其中，l为每个单元均匀圆阵Adcock天线的圆阵半径，λ为外辐射源的调频广播信号的波长；

(5.2)对子步骤(5.1)中建立的目标方位角的优化模型进行求解，得出目标方位角θ_s。

6.如权利要求5所述的基于调频广播信号的外辐射源雷达修正比相测角方法，其特征在于，在子步骤(5.2)中，对子步骤(5.1)中建立的目标方位角的优化模型进行求解的过程包括以下子步骤：

a)根据公式 $\left[\begin{array}{c}{T}^{*}\\ {\mathrm{\&xi;}}^{*}\end{array}\right]={\left(A^H\left(\mathrm{\&theta;}\right)A\left(\mathrm{\&theta;}\right)\right)}^{-1}{A}^H\left(\mathrm{\&theta;}\right)R,$ 计算数值变量T与ξ的最小二乘解，其中，上标*表示共轭，上标H表示矩阵的共轭转置，上标-1表示矩阵的逆，A(θ)＝[a(θ)-b]；将求解出的T与ξ的最小二乘解代入ψ(θ,T,ξ)中，得出函数ψ(θ)；

b)设定角度初始值θ₁，设定允许误差

c)设置参数H₁＝1，计算ψ(θ)在角度初始值θ₁处的梯度g₁，设置整数变量m，m＝1,2...，当m＝1时，执行子步骤d)；

d)计算搜索方向d_m，d_m＝-H_mg_m；

e)求解迭代步长λ_m，使其满足 $\mathrm{\&psi;}({\mathrm{\&theta;}}_m+{\mathrm{\&lambda;}}_m{d}_m)=\underset{{\mathrm{\&lambda;}}^{\&prime;}\&GreaterEqual;0}{\min }\mathrm{\&psi;}({\mathrm{\&theta;}}_m+{\mathrm{\&lambda;}}^{\&prime;}{d}_m),$ 其中，λ'为大于等于0的变量，表示：当λ'在[0,+∞)变化时ψ(θ_m+λ'd_m)能取到的最大值；然后，令θ_m+1＝θ_m+λ_md_m；

f)判断与允许误差的大小关系，表示ψ(θ)在θ_m+1处的梯度，表示的绝对值；若则得出θ_s＝θ_m+1；否则进行子步骤g)；

g)令p_m＝θ_m+1-θ_m，q_m＝g_m+1-g_m，然后计算H_m+1，

$H_{m+1}=H_m+\frac{p_m{p}_m}{p_m{q}_m}-\frac{H_m{q}_m{q}_m{H}_m}{q_m{H}_m{q}_m}$

然后，令m的值自增1，返回子步骤d)。