CN104515971B - 宽带多目标机载单站无源定位方法 - Google Patents

Abstract

本发明公开了一种宽带多目标机载单站无源定位方法，主要解决现有技术不能对多宽带辐射源目标进行单站无源定位的问题。其实现步骤是：(1)在接收端对远场宽带辐射源信号进行采样，得到观测信号；(2)通过分数阶傅里叶变换FRFT将阵元接收的宽带观测信号窄带化；(3)构建第k个信号在分数阶傅里叶变换FRFT域中的信号向量I^(k)；(4)对第k个信号利用信号向量I^(k)求时间平均自相关矩阵，并做特征值分解；(5)构建第k个信号距离形式的导向矩阵；(6)在接收端进行L次采样，重复步骤(2)至(5)，对目标函数搜索极大值，极大值点处对应的(x_k,y_k,z_k)为目标k的坐标。本发明能对多个宽带辐射源目标进行单站无源定位，可用于无人、载人飞行平台探测侦察。

Description

宽带多目标机载单站无源定位方法

技术领域

本发明属于雷达技术领域，特别涉及一种宽带多目标机载单站无源定位方法，可用于电子侦察监视、远距离预警探测领域多辐射源目标的无源定位。

背景技术

有源装备如雷达、声纳等，主动发射电磁波信号，接收来自目标反射的回波信号，此回波信号含有目标位置信息，进一步解算即可实现目标定位。这种通过装备自身发射电磁信号实现目标定位的方法称为有源定位。与有源定位相比，不主动发射信号而被动接收单程的目标辐射源信号的无源定位具有作用距离远，隐蔽性强的优势，在反辐射打击、低空突防等威胁力越来越强的现代战场环境中，无源定位系统具有更强的生存能力。

通过单个观测站截获非合作辐射源信号、获取目标信息，并以一定精度给出目标空间坐标的单站无源定位，具有较好的抗反辐射打击、抗低空突防的能力。与分布式多站无源定位相比，单站无源定位系统的战术机动性和对诸如海基、空基等定位场景的适应性更强，克服了多站系统在复杂战场环境下的几何布站、站间通信等方面的棘手问题，成为无源雷达技术领域中的研究重点。

传统的单站无源定位法中，只测向BOL定位法对方向测量误差敏感，对观测站机动性要求高，定位实时性差、定位精度低。目前研究较为集中的基于质点运动学的测空、时、频域参量变化率信息，即方位角速度、俯仰角速度、相位变化率、脉冲重复周期变化率和多普勒频率变化率的单站无源定位法，定位速度快、定位精度高，但对参量变化率的测量精度要求较高；此外，解定位模糊也是该类定位法中一个难点。例如，张敏,郭福成,周一宇，“基于单个长基线干涉仪的运动单站直接定位”，航空学报，2013,34(2),pp.378-386，提出通过观测站匀速运动给出多个定位圆以解算目标位置，实现无模糊定位，但若要提高收敛速度除强机动性以外，对观测站运动姿态也提出一定要求，以使定位线簇之间有较明显差异，否则将形成多条模糊定位区域；还有一部分学者利用滤波算法和代价函数解定位模糊，例如，郭福成,贾兴江,皇甫堪，“仅用相位差变化率的机载单站无源定位方法及其误差分析”，航空学报，2009,30(6),pp.1090-1095，提出采用最小二乘、扩展卡尔曼(EKF)等非线性跟踪滤波方法求解目标位置，但该算法求解过程中需要给定初始值，初始值的选取及其估计精度直接影响定位结果。以上传统定位法需要先单个或多个参量估计，继而求解定位代价函数，定位参数与目标辐射源频率大小有关，对天线的方向性及其变化率有较高要求，同时在观测站的机动性和运动姿态等方面也存在很多限制；此外，针对的辐射源信号均为窄带，对宽带场景不适用。

发明内容

本发明的目的在于克服上述已有技术的不足，提出宽带多目标机载单站无源定位方法，以在无需参数测量的条件下，从宽带信号数据中直接解析出目标位置，并在低信噪比下提高定位性能。

本发明的技术方案是：通过在接收端对远场宽带辐射源信号采样，得到观测信号；通过分数阶傅里叶变换FRFT将阵元接收的宽带观测信号窄带化；利用窄带化后的信号构建信号向量；利用信号向量求时间平均自相关矩阵并做特征值分解；构建目标函数中的距离形式的导向矩阵；L次测量在目标位置处累积去模糊，对目标函数搜索极大值，以实现目标定位。其实现步骤如下:

(1)在接收端对远场宽带辐射源信号进行采样，得到观测信号；

(2)通过分数阶傅里叶变换FRFT将阵元接收的宽带观测信号窄带化；

(3)构建第k个信号在分数阶傅里叶变换FRFT域中的信号向量I^(k)，其中，k＝1,2,...,Q，Q为辐射源数；

(4)对第k个信号利用信号向量I^(k)求时间平均自相关矩阵，并对时间平均自相关矩阵做特征值分解；

(5)构建目标函数中第k个信号距离形式的导向矩阵A_k,l，其中，l为采样次数序号，l＝1,2,...,L，L为总的采样次数；

(6)在接收端进行L次采样，重复步骤(2)至(5)，搜索目标函数的极大值，极大值点处对应的(x_k,y_k,z_k)为目标k的空间坐标，其中，v表示噪声，为第l次采样时得到的噪声子空间矩阵，H表示取共轭转置。

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

第一，本发明利用观测信号数据直接解析出目标位置，克服了现有技术中定位参数估计困难的缺点，使得本发明对不同定位场景和搭载平台都具有较好的适应性。

第二，本发明利用分数阶傅里叶变换的线性性质及其对宽带信号的能量聚焦性，在无交叉项影响下将宽带信号窄带化，实现了多宽带辐射源信号的检测与定位问题，克服了现有定位技术不适用于宽带定位场景的缺点。

第三，本发明中目标函数的极大值搜索在距离维进行，从而可以在低信噪比和较小采样快拍数时也能利用多维数组的极值搜索算法较为准确地搜索到各极大值，取得较好的定位性能。

附图说明

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

图2为本发明的宽带多目标机载单站无源定位原理图；

图3为本发明中均匀圆阵天线图；

图4为宽带辐射源定位场景图；

图5为信噪比5dB时0号阵元上观测信号的离散分数阶傅里叶变换DFRFT图；

图6为信号比5dB时对两个目标A、B的定位结果图；

图7为在不同信噪比下两个目标A、B的X、Y坐标均方根误差曲线图。

具体实施方式

以下参照附图对本发明作进一步详细的描述。

参照图1，本发明实现步骤如下：

步骤1：在接收端对远场宽带辐射源信号进行采样，得到观测信号。

(1.1)观测站沿设定航迹飞行，其航迹参数可通过自身设备获得，设接收端的天线为均匀圆阵天线，如图3所示，该天线的半径为d，阵元数为M，阵元编号m＝0,1,...,M-1，阵元参考点为圆心；

(1.2)利用上述接收端的天线接收远场宽带辐射源信号，对接收到的信号做等间隔采样，得到观测信号其中，s_k为第k个远场宽带辐射源信号，k＝1,2,...,Q，n＝1,2,...,N，N为采样快拍数，V^m(n)为加性高斯白噪声，为第k个信号到达阵元m时相对于参考点的时延。

步骤2：通过分数阶傅里叶变换FRFT将阵元接收的宽带观测信号窄带化。

(2.1)将第0号阵元上接收的第k个信号s_k(n)表示为：

其中，f_s为采样率，f_k和μ_k分别为信号s_k(n)的初始频率和调频斜率；

(2.2)对上述信号s_k(n)做离散分数阶傅里叶变换，得到变换后的信号

其中，K_α(n,ε)为离散分数阶傅里叶变换的核函数，p为离散分数阶傅里叶变换的阶数，α＝pπ/2为旋转因子，ε为信号s_k(n)在分数阶傅里叶变换域中的频点；

(2.3)对变换后的信号进行二维搜索，得到极大值为：

其中，最佳旋转因子

(2.4)根据最佳旋转因子α_k对第m个阵元接收到的第k个信号做旋转因子为α_k的离散分数阶傅里叶变换，得到变换量W^(m,k)(α_k,ε)：

其中，为第m个阵元接收到的第k个信号，{F_α[·]}(ε)表示做离散分数阶傅里叶变换；

(2.5)对变换量W^(m,k)(α_k,ε)进行搜索，得到变换量的极大值W^(m,k)(α_k,ε_k,m)：

其中，

步骤3：根据步骤2的结果构建第k个信号在分数阶傅里叶变换FRFT域中的信号向量I^(k)，按如下公式构建：

I^(k)＝[W^(0,k)(α_k,ε_k,0),W^(1,k)(α_k,ε_k,1),...,W^(m,k)(α_k,ε_k,m),...,W^(M-1,k)(α_k,ε_k,M-1)]^T，

其中，W^(m,k)(α_k,ε_k,m)为向量I^(k)的元素，m＝0,1,...,M-1，T表示取转置。

步骤4：对第k个信号利用信号向量I^(k)求时间平均自相关矩阵，并对时间平均自相关矩阵做特征值分解。

(4.1)计算第k个信号的时间平均自相关矩阵R_k：

R_k＝E{I^(k)[I^(k)]^H}，

其中，E{·}表示取均值；

(4.2)对时间平均自相关矩阵R_k做特征值分解：

其中，v表示噪声，s表示信号，G^s、G^v分别表示信号子空间矩阵和噪声子空间矩阵，Λ为特征值组成的特征矩阵。

步骤5：构建目标函数中第k个信号距离形式的导向矩阵A_k,l。

(5.1)根据图2中给出的观测站与远场辐射源目标的几何关系，计算接收端的单位距离矢量ρ_k,l：

其中，△r_k,l＝Γ_k-O_l＝[△x_k,l,△y_k,l,△z_k,l]为第l次采样时观测站与目标k之间的距离矢量，Γ_k＝(x_k,y_k,z_k)为目标k的坐标，O_l＝(x_obs,y_obs,z_obs)为观测站的坐标，obs表示观测站，△x_k,l、△y_k,l、△z_k,l分别为目标k的坐标和观测站的坐标中对应分量的差值，λ_k＝c/f_k，c为电波传输速度，|·|表示取模；

(5.2)利用单位距离矢量ρ_k,l构建目标函数中第k个信号距离形式的导向矩阵A_k,l：

其中，为导向矩阵A_k,l的元素，m＝0,1,…,M-1，为第l次采样时第k个信号到达阵元m时相对于参考点的时延，η_m＝[dcos(2πm/M),dsin(2πm/M),0]^T为均匀圆阵天线中第m号阵元的坐标，为第k个目标的来波入射方向的单位向量，θ_k和分别为第k个目标的俯仰角和方位角。

步骤6：在接收端进行L次采样，重复步骤2至5，对目标函数搜索极大值，极大值点处对应的(x_k,y_k,z_k)为目标k的空间坐标。

所述目标函数表示如下：

其中，v表示噪声，为第l次采样时得到的噪声子空间矩阵。

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

1)仿真场景与参数

为验证本发明的有效性和正确性，给出了宽带辐射源定位场景如图4所示，定位场景中的观测站初始坐标为(0,0,10000)，两个远场辐射源目标A、B的坐标分别为(10000,2000,0)、(8000,-2000,0)，坐标单位均为米。观测站速度300m/s，做‘S’形运动，采样次数取10，采样率取100MHz，采样快拍数取512。

所有仿真实验均在Windows 7操作系统下采用Matlab 2010a软件实现。

设目标A的初始频率为10MHz，目标A的调频斜率为8MHz/μs；设目标B的初始频率为30MHz，目标B的调频斜率为-8MHz/μs；均匀圆阵天线的阵元数为10，如图3所示。

为了评价定位方法的性能，采用如下均方根误差RMSE作为评价因子：

其中，为辐射源目标的坐标估计值，X、Y为辐射源目标的坐标真实值，E{·}表示取均值。

2)仿真内容

通过本发明方法对场景中的两个目标进行定位。

仿真1：在信噪比5dB时，用本发明对0号阵元上的观测信号做离散分数阶傅里叶变换DFRFT，得到两个目标A、B的变换峰，峰值点坐标分别为(p_A,ε_A)和(p_B,ε_B)，结果如图5所示。

根据变换峰的坐标参数得到目标函数中的导向矩阵和噪声子空间矩阵，并对目标函数搜索其极大值实现对两个目标A、B的定位，结果如图6所示。

由图5的仿真计算结果可知，得到的峰值点坐标(p_A,ε_A)和(p_B,ε_B)均正确。

由图6的定位结果可知，在信噪比5dB时，用本发明实现了两个目标A、B的定位，其中，目标A的X、Y坐标估计值分别为9500m和2300m，如图6(a)所示；目标B的X、Y坐标估计值分别为8000m和-2200m，如图6(b)所示。

仿真2：在信噪比变化为0-15dB的范围，做100次MonteCarlo仿真实验，得到两个目标A、B的X坐标均方根误差随信噪比的变化曲线，如图7(a)所示；目标A、B的Y坐标均方根误差随信噪比的变化曲线，如图7(b)所示。

由图7可知，在信噪比变化为0-15dB的范围，本发明均能实现两个宽带辐射源目标的检测与定位，且在信噪比为0dB时仍具有较好的定位性能。

Claims (3)

1.一种宽带多目标机载单站无源定位方法，包括如下步骤：

(1)在接收端对远场宽带辐射源信号进行采样，得到观测信号；

(2)通过分数阶傅里叶变换FRFT将阵元接收的宽带观测信号窄带化:

(2.1)设接收端的天线为均匀圆阵天线，其半径为d，阵元数为M，阵元编号依次为0,1,…,M-1，阵元参考点为圆心，将第0号阵元上接收的第k个信号s_k(n)表示为：

$s_k\left(n\right)=e^{j\mathrm{\π}\[2f_{k}n/f_s+{\mathrm{\μ}}_k{(n/f_s)}^2\]},$

其中，f_s为采样率，N为采样快拍数，n＝1,2,…,N，f_k和μ_k分别为信号s_k(n)的初始频率和调频斜率；

(2.2)对上述s_k(n)做离散分数阶傅里叶变换，得到变换后的信号

$\begin{array}{c}{X}_p^{\left(k\right)}(\mathrm{\α},\mathrm{\ϵ})=\underset{n=-(N-1)/2}{\overset{(N-1)/2}{\mathrm{\Σ}}}{K}_{\mathrm{\α}}(n,\mathrm{\ϵ})s_k\left(n\right)\\ =\frac{\sqrt{1-j\cot \mathrm{\α}}}{\sqrt{N}}{e}^{{\mathrm{j\π\ϵ}}^2\cot \mathrm{\α}/N}\underset{n=-(N-1)/2}{\overset{(N-1)/2}{\mathrm{\Σ}}}\[e^{j\mathrm{\π}n(-2\mathrm{\ϵ}\csc \mathrm{\α}/N+2f_k/f_s}{e}^{{\mathrm{j\πn}}^2(\cot \mathrm{\α}/N+{\mathrm{\μ}}_k/f_s^2)}\]\end{array},$

其中，K_α(n,ε)为离散分数阶傅里叶变换的核函数，p为离散分数阶傅里叶变换的阶数，α＝pπ/2为旋转因子，ε为信号s_k(n)在分数阶傅里叶变换域中的频点；

(2.3)对变换后的信号进行二维搜索，得到极大值为：

$X_p^{\left(k\right)}({\mathrm{\α}}_k,{\mathrm{\ϵ}}_k)=\sqrt{N(1-{\mathrm{jcot\α}}_k)}{e}^{{\mathrm{j\π\ϵ}}_k^2{\mathrm{cot\α}}_k/N},$

其中，最佳旋转因子ε_k＝f_kNsinα_k/f_s；

(2.4)根据最佳旋转因子α_k对第m个阵元接收到的第k个信号做旋转因子为α_k的离散分数阶傅里叶变换，得到变换量W^(m,k)(α_k,ε)：

$W^{(m,k)}({\mathrm{\α}}_k,\mathrm{\ϵ})=\{F_{{\mathrm{\α}}_k}\[s_k(n-{\mathrm{\τ}}_k^m)\]\}\left(\mathrm{\ϵ}\right),$

其中，m＝0,1,…,M-1，为第m个阵元接收到的第k个信号，为信号到达阵元m时相对于参考点的时延，{F_α[·]}(ε)表示做离散分数阶傅里叶变换；

(2.5)对变换量W^(m，k)(α_k,ε)进行搜索，得到变换量的极大值W^(m,k)(α_k，ε_k，m)：

$\begin{array}{c}{W}^{(m,k)}{({\mathrm{\α}}_k,{\mathrm{\ϵ}}_{k,m})}^{j\mathrm{\π}\[-2f_k{\mathrm{\τ}}_k^m+{\mathrm{\μ}}_k{\left({\mathrm{\τ}}_k^m\right)}^2\]}\sqrt{N(1-{\mathrm{jcot\α}}_k)}{e}^{j\mathrm{\π}{\left({\mathrm{\ϵ}}_{k,m}\right)}^2{\mathrm{cot\α}}_k/N}\\ =U_k\left({\mathrm{\τ}}_k^m\right)X_p^{\left(k\right)}({\mathrm{\α}}_k,{\mathrm{\ϵ}}_k)\end{array},$

其中，

(3)构建第k个信号在分数阶傅里叶变换FRFT域中的信号向量I^(k)，其中，k＝1,2,…,Q，Q为辐射源数；

(4)对第k个信号利用信号向量I^(k)求时间平均自相关矩阵，并对时间平均自相关矩阵做特征值分解；

(5)构建目标函数中第k个信号距离形式的导向矩阵A_k,l，其中，l为采样次数序号，l＝1,2,…,L，L为总的采样次数；

(6)在接收端进行L次采样，重复步骤(2)至(5)，搜索目标函数的极大值，极大值点处对应的(x_k,y_k,z_k)为目标k的空间坐标，其中，v表示噪声，为第l次采样时得到的噪声子空间矩阵，H表示取共轭转置。

2.根据权利要求1所述的方法，其中步骤(3)所述的构建第k个信号在分数阶傅里叶变换FRFT域中的信号向量I^(k)，按如下公式构建：

I^(k)＝[W^(0,k)(α_k，ε_k,0),W^(1,k)(α_k，ε_k,1),…,W^(m,k)(α_k，ε_k,m),…,W^(M-1,k)(α_k，ε_k,M-1)]^T，

其中，W^(m,k)(α_k，ε_k,m)为向量I^(k)的元素，m＝0,1,…,M-1，T表示取转置。

3.根据权利要求1所述的方法，其中步骤(5)所述的构建目标函数中第k个信号距离形式的导向矩阵A_k,l，按如下步骤进行：

(5.1)计算接收端的单位距离矢量ρ_k,l：

${\mathrm{\ρ}}_{k,l}=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_k}\frac{{\mathrm{\Δr}}_{k,l}}{\left|{\mathrm{\Δr}}_{k,l}\right|},$

其中，Δr_k,l＝Γ_k-O_l＝[Δx_k,l,Δy_k,l,Δz_k,l]为第l次采样时观测站与目标k之间的距离矢量，Γ_k＝(x_k,y_k,z_k)为目标k的坐标，O_l＝(x_obs,y_obs,z_obs)为观测站的坐标，obs表示观测站，Δx_k,l、Δy_k,l、Δz_k,l分别为目标k的坐标和观测站的坐标中对应分量的差值，λ_k＝c/f_k，c为电波传输速度，|·|表示取模；

(5.2)利用单位距离矢量ρ_k,l构建目标函数中第k个信号距离形式的导向矩阵A_k,l：

$A_{k,l}=\[U_k\left({\mathrm{\τ}}_{k,l}^0\right),U_k\left({\mathrm{\τ}}_{k,l}^1\right),...,U_k\left({\mathrm{\τ}}_{k,l}^m\right),...,U_k\left({\mathrm{\τ}}_{k,l}^{M-1}\right)\],$

其中，为导向矩阵A_k,l的元素，m＝0,1,…,M-1，为第l次采样时第k个信号到达阵元m时相对于参考点的时延，η_m＝[d cos(2πm/M),dsin(2πm/M),0]^T为均匀圆阵天线中第m号阵元的坐标，为第k个目标的来波入射方向的单位向量，θ_k和分别为第k个目标的俯仰角和方位角。