在现代电子战中对未知目标的无源定位是当前研究的主要问题之一，在无源定位系统中观测站本身不辐射电磁波，基本为完全被动工作，因此，具有良好的隐蔽性、抗干扰能力强的优点，备受各国的关注和重视^[1-2]。无源定位通常首先通过现有的测量技术获取辐射源信号的参数；然后，利用不同的参数采用相应定位方法^[3-4]；最后，建立观测模型，根据观测模型选择正确的定位算法。

多站时差定位技术是通过多个观测站对未知辐射源信号进行信号信息参数的收集，然后通过相应的处理获取TDOA进行时差定位^[4-7]。时差定位的定位精度比较高，但是，也容易出现时差模糊的情况，对于高重频的目标信号定位难度比较大，并且对于目标速度无法确定。在时差定位的基础上加入多普勒频差信息能提高定位精度，消除时差模糊的问题，而且对目标的速度也能确定^[7-8]。

利用标准粒子群算法对辐射源定位容易出现局部最优解，导致定位失败^[8-12]。针对这种情况本文利用自然选择和遗传学原理进行选择淘汰和交叉、变异，对下一代粒子进行变化来达到种群的多样性使定位解算能够快速地收敛以及达到全局最优。

1 时差频差联合定位原理

本文是在三维空间里进行定位，取地固直角坐标系为基准坐标系，如图 1所示。设辐射源T的位置矢量为 U=[x, y, z]，速度矢量为$\dot{\boldsymbol{U}}=[\dot{x}, \dot{y}, \dot{z}]$，观测站S_i坐标矢量，速度矢量${\mathit{\boldsymbol{\dot X}}_i} = \left[ {{{\dot x}_i},{{\dot y}_i},{{\dot z}_i}} \right]$], i=(1, 2, …, M)。要对运动目标位置和速度进行定位，至少需要M≥4个观测站同时产生不少于3对到达时间差(time difference of arrival, TDOA)、到达频差(frequence difference of arrival, FODA)测量值。

不失一般性，本文不考虑4个观测站在同一直线和同一平面的情况。假设S₁作为主站，辐射源到各个观测站的距离为:

$ \begin{array}{*{20}{c}} {{R_i} = \left\| {\mathit{\boldsymbol{U}} - {\mathit{\boldsymbol{X}}_i}} \right\| = \sqrt {{{\left( {{x_i} - x} \right)}^2} - {{\left( {{y_i} - y} \right)}^2} - {{\left( {{z_i} - z} \right)}^2}} }\\ {(i = 1,2,3,4)} \end{array} $

(1)

式中‖·‖为2的范数。

辐射源到观测站1和到观测站i的距离差为

$ {R_{i1}} = {R_i} - {R_1} = c{t_{i1}}(i = 2,3, \cdots ,M) $

(2)

式中：c为光速；t_i1为到达时间差。

对式(2)时间求微分可得到一组多普勒频差观测方程：

$ {{\dot R}_{i1}} = {{\dot R}_i} - {{\dot R}_1} = \frac{{c{f_{i1}}}}{{{f_c}}}(i = 2,3, \cdots ,M) $

(3)

式中：f_i1为多普勒频差；f_c为载波频率；Ṙ_i为R_i的变化率。对式(1)时间求微分可得：

$ {\dot R_i} = \frac{{{{\left( {\mathit{\boldsymbol{\dot U}} - {{\mathit{\boldsymbol{\dot X}}}_i}} \right)}^{\rm{T}}}\left( {\mathit{\boldsymbol{U}} - {\mathit{\boldsymbol{X}}_i}} \right)}}{{{R_i}}} $

(4)

在实际工程中，观测量TDOA=[t₂₁, t₃₁, …, t_M1]、FDOA=[f₂₁, f₃₁, …, f_M1]，一般都存在加性噪声，所以得到的时差与频差观测方程为

$ \left\{ {\begin{array}{*{20}{l}} {R_{i1}^o = c{t_{i1}} = {R_{i1}} + {n_{i1}}}\\ {\dot R_{i1}^o = \frac{{c{f_{i1}}}}{{{f_c}}} = {{\dot R}_{i1}} + {{\dot n}_{i1}}} \end{array}} \right. $

(5)

式中：R_i1^o和Ṙ_i1^o分别为含噪声的距离差和多普勒频差；n_i1和ṅ_i1分别为对应的测量误差。

可将时差与频差的观测方程改为向量形式为：

$ \mathit{\boldsymbol{\varepsilon }} = \mathit{\boldsymbol{h}}(\mathit{\boldsymbol{\theta }}) + \mathit{\boldsymbol{n}} $

(6)

式中：ε= [R_i1^o, Ṙ_i1^o]^T为时差和频差的观测值信息；$\mathit{\boldsymbol{\theta }} = {[U,\dot U]^{\rm{T}}}$辐射源位置和速度信息；h(θ)=[R_i1, Ṙ_i1]^T可看作要求解的未知量θ与已知量R_i1、Ṙ_i1之间存在非线性的函数关系；n 为均值为零、协方差矩阵为Q高斯噪声。

2 优化PSO算法的时差频差联合定位算法

由于标准粒子群算法存在早期收敛，导致收敛精度不高，并且很容易有局部最优的问题。其中遗传算法中是利用选择、交叉和变异算子进行搜索，全局搜索的能力比较强，局部搜索的能力比较弱，因此，结合PSO算法, 利用粒子群算法进行局部搜索，利用遗传算法进行全局搜索，从而得到辐射源的全局最优解。因此，本节对PSO进行了改进来提高其收敛速度和精度，并且使其避开局部最优的问题，达到全局最优的目的。

在进行定位求解的问题时，就是考虑最小化的问题，首先，设置空间搜索维数为D维，群体规模为N，粒子i第t时刻的位置和速度分别为：

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{X}}_i}(t) = \left( {{X_{i1}}(t),{X_{i2}}(t), \cdots ,{X_{iD}}(t)} \right)}\\ {{\mathit{\boldsymbol{V}}_i}(t) = \left( {{V_{i1}}(t),{V_{i2}}(t), \cdots ,{V_{iD}}(t)} \right)} \end{array}} \right. $

(7)

第i个粒子的自身最优位置P_best为P_i(t)=(P_i1(t), P_i2(t), …, P_iD(t))；群体的最优位置G_best为G_i(t)=(G_i1(t), G_i2(t), …, G_iD(t))。粒子的速度与位置信息按照式(8)进行更新：

$ \left\{ \begin{array}{l} {V_{id}}(t + 1) = w(t){V_{id}}(t) + {c_1}r{1_{id}}(t)\left( {{P_{id}}(t)} \right.\\ \left. {\;\;\;\;\;\; - {X_{id}}(t)} \right) + {c_2}{r_2}(t)\left( {{G_{id}}(t) - {X_{id}}(t)} \right)\\ {X_{id}}(t + 1) = {X_{id}}(t) + {V_{id}}(t + 1) \end{array} \right. $

(8)

式中:r1和r2为[0, 1]的相互独立，且均匀分布随机数; c₁、c₂学习因子，值视情况而定。

w为线性惯性权重值，其更新策略采用收敛速度更快的方式：

$ w(t) = {w_{{\rm{end}}}} + \left( {{w_{{\rm{init}}}} - {w_{{\rm{end}}}}} \right)\left( {{t_{{\rm{max}}}} - t} \right)/{t_{{\rm{max}}}} $

(9)

式中：w_end代表最大权重因子；w_init代表最小权重因子；t_max为最大迭代次数；w_end、w_init分别为0.8、0.2。在带有惯性权重系数的PSO算法中，w值随着迭代次数的增加而减少。在迭代前期，较小的w值利于进化前期加强探测；在迭代后期，较小的w值利于对当前优秀解的开发，以便得到更为精确的解。

在得出新的位置和速度信息后，将速度和位置限制为一定范围内：其中位置X_ij(t)∈[-X_max, X_min]，速度V_ij(t)∈[-V_max, V_min]。当有粒子超过这个范围时，重新更新该粒子对应的信息，更新方式为在这个范围内随机更新。

粒子i更新位置后，需要更新自身的P_best值，其更新方程为：

$ {P_i}(t + 1) = \left\{ {\begin{array}{*{20}{l}} {{X_i}(t + 1),f\left( {{X_i}(t + 1)} \right) < f\left( {{P_i}(t)} \right)}\\ {{P_i}(t),a} \end{array}} \right. $

(10)

式中：f为适应度函数，对于最小化问题，其值越小、则对应的解越优秀。

当群体中所有粒子的P_best更新后，将其最优粒子G_best进行更新，其更新方程为：

$ G(t + 1) = \left\{ {\begin{array}{*{20}{l}} {{P_k}(t + 1),f\left( {{P_k}(t + 1)} \right) < f(G(t))}\\ {G(t),其他} \end{array}} \right. $

(11)

式中$k = \arg {\min\limits_{1 \le i \le N}}\left\{ {f\left({{P_i}(t)} \right)} \right\}$，即适应度值最小的粒子。

将更新后的种群按照适应值排序，用群体最好的一半粒子的位置和速度替代最差的一半粒子，并且保持P_best和G_best不变。

然后，根据杂交概率选取一定数量的粒子放入杂交池中，在杂交池中的粒子两两进行杂交产生相同数目的子代粒子，并保持P_best和G_best不变，其中子代粒子的位置和速度更新公式为

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{child}}(x) = p \cdot {\rm{paren}}{{\rm{t}}_1}(x) + (1 - p) \cdot {\rm{paren}}{{\rm{t}}_2}(x)}\\ {{\rm{child}}\;(v) = \frac{{{\rm{paren}}{{\rm{t}}_1}\;(v)\;+\;{\rm{paren}}{{\rm{t}}_2}\;(v)}}{{|{\rm{paren}}{{\rm{t}}_1}\;(v)\;+\;{\rm{paren}}{{\rm{t}}_2}\;(v)|}}\;\;\;\;\;\;|{\rm{paren}}{{\rm{t}}_2}\;(v)|} \end{array}} \right. $

(12)

式中：child为子代粒子的信息；parent为父代粒子的信息。最后输出结果及辐射源的位置和速度信息。

在确定PSO算法之后，要确定适应度函数，适应度函数的推导如下。

根据时差方程式(2)可以得到：

$ {R_i} = {R_{i1}} + {R_1} $

(13)

将式(13)两边进行平方移项后得出新的一组时差方程为

$ \begin{array}{*{20}{c}} {R_{i1}^2 + 2{R_{i1}}{R_1} = R_i^2 - R_1^2 = }\\ {{{\left( {{\mathit{\boldsymbol{X}}_i} - \mathit{\boldsymbol{U}}} \right)}^{\rm{T}}}\left( {{\mathit{\boldsymbol{X}}_i} - \mathit{\boldsymbol{U}}} \right) - {{\left( {{\mathit{\boldsymbol{X}}_1} - \mathit{\boldsymbol{U}}} \right)}^{\rm{T}}}\left( {{\mathit{\boldsymbol{X}}_1} - \mathit{\boldsymbol{U}}} \right) = }\\ {\mathit{\boldsymbol{X}}_i^{\rm{T}}{\mathit{\boldsymbol{X}}_i} - \mathit{\boldsymbol{X}}_1^{\rm{T}}{\mathit{\boldsymbol{X}}_1} - 2{{\left( {{\mathit{\boldsymbol{X}}_i} - {\mathit{\boldsymbol{X}}_1}} \right)}^{\rm{T}}}\mathit{\boldsymbol{U}},}\\ {i = 2,3, \cdots ,\mathit{\boldsymbol{M}}} \end{array} $

(14)

由于式(5)中含有噪声的时差方程，结合式(14)可得：

$ {\varepsilon _{ti}} = R_{i1}^2 + 2{R_{i1}}{R_1} - \mathit{\boldsymbol{X}}_i^{\rm{T}}{\mathit{\boldsymbol{X}}_i} + \mathit{\boldsymbol{X}}_1^{\rm{T}}{\mathit{\boldsymbol{X}}_1} + 2{\left( {{\mathit{\boldsymbol{X}}_i} - {\mathit{\boldsymbol{X}}_1}} \right)^{\rm{T}}}\mathit{\boldsymbol{X}} $

(15)

则可构造时差方程的矩阵为:

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\varepsilon }}_t} = {\mathit{\boldsymbol{h}}_t} - {\mathit{\boldsymbol{G}}_t}\mathit{\boldsymbol{\theta }}}\\ {{\mathit{\boldsymbol{\varepsilon }}_t} = \left[ {{\mathit{\boldsymbol{\varepsilon }}_{t1}},{\mathit{\boldsymbol{\varepsilon }}_{t2}}, \cdots ,{\mathit{\boldsymbol{\varepsilon }}_{tM}}} \right]} \end{array} $

(16)

$ {\mathit{\boldsymbol{h}}_t} = \left[ {\begin{array}{*{20}{c}} {R_{21}^2 - \mathit{\boldsymbol{X}}_2^{\rm{T}}{\mathit{\boldsymbol{X}}_2} + \mathit{\boldsymbol{X}}_1^{\rm{T}}{\mathit{\boldsymbol{X}}_1}}\\ {\mathit{\boldsymbol{R}}_{31}^2 - \mathit{\boldsymbol{X}}_3^{\rm{T}}{\mathit{\boldsymbol{X}}_3} + \mathit{\boldsymbol{X}}_1^{\rm{T}}{\mathit{\boldsymbol{X}}_1}}\\ \vdots \\ {R_{M1}^2 - \mathit{\boldsymbol{X}}_M^{\rm{T}}{\mathit{\boldsymbol{X}}_M} + \mathit{\boldsymbol{X}}_1^{\rm{T}}{\mathit{\boldsymbol{X}}_1}} \end{array}} \right] $

$ {\mathit{\boldsymbol{G}}_t} = - 2\left[ {\begin{array}{*{20}{c}} {{{\left( {{\mathit{\boldsymbol{X}}_2} - {\mathit{\boldsymbol{X}}_1}} \right)}^{\rm{T}}}}&{{R_{21}}}&{{{\bf{0}}^{\rm{T}}}}&0\\ {{{\left( {{\mathit{\boldsymbol{X}}_3} - {\mathit{\boldsymbol{X}}_1}} \right)}^{\rm{T}}}}&{{R_{31}}}&{{{\bf{0}}^{\rm{T}}}}&0\\ {}& \vdots &{}&{}\\ {{{\left( {{\mathit{\boldsymbol{X}}_M} - {\mathit{\boldsymbol{X}}_1}} \right)}^{\rm{T}}}}&{{R_{{\rm{M1}}}}}&{{{\bf{0}}^{\rm{T}}}}&0 \end{array}} \right] $

$\mathit{\boldsymbol{\theta }} = \left[{\mathit{\boldsymbol{U}}, {R_1}, \mathit{\boldsymbol{\dot U}}, {{\dot R}_1}} \right]$，其中0为3行1列的零矩阵。

式(14)是一个含有辐射源U位置信息的非线性方程，对于无源定位中只含时差信息观测参数只能对辐射源位置进行解算，通过引入FDOA信息与TDOA结合，不仅能解出目标的瞬时速度信息，还能获得比单纯TDOA定位更高的位置精度。因此，对式(14)进行求导可得：

$ \begin{array}{*{20}{c}} {{R_{i1}}{{\dot R}_{i1}} + {{\dot R}_{i1}}{R_1} + {R_{i1}}{{\dot R}_1} = }\\ {\mathit{\boldsymbol{X}}_i^{\rm{T}}{\mathit{\boldsymbol{X}}_i} - \mathit{\boldsymbol{\dot X}}_1^{\rm{T}}{\mathit{\boldsymbol{X}}_1} - {{\left( {\mathit{\boldsymbol{\dot X}}_1^{\rm{T}} - {\mathit{\boldsymbol{X}}_1}} \right)}^{\rm{T}}}\mathit{\boldsymbol{U}} - {{\left( {\mathit{\boldsymbol{X}}_i^{\rm{T}} - {\mathit{\boldsymbol{X}}_1}} \right)}^{\rm{T}}}\mathit{\boldsymbol{\dot U}}} \end{array} $

结合式(5)中频差的测量方程可得时差方程的矩阵形式为：

$ \begin{array}{l} \;\;\;\;\;\;{\mathit{\boldsymbol{\varepsilon }}_f} = {\mathit{\boldsymbol{h}}_f} - {\mathit{\boldsymbol{G}}_f}\mathit{\boldsymbol{\theta }}\\ {\mathit{\boldsymbol{\varepsilon }}_f} = \left[ {{\mathit{\boldsymbol{\varepsilon }}_{f1}},{\varepsilon _{f2}}, \cdots ,{\varepsilon _{fM}}} \right] \end{array} $

(17)

$\boldsymbol{h}_{f}=2\left[\begin{array}{c}{R_{21} \dot{R}_{21}-\dot{\boldsymbol{X}}_{2}^{\mathrm{T}} \boldsymbol{X}_{2}+\dot{\boldsymbol{X}}_{1}^{\mathrm{T}} \boldsymbol{X}_{1}} \\ {R_{31} \dot{R}_{31}-\dot{\boldsymbol{X}}_{3}^{\mathrm{T}} \boldsymbol{X}_{3}+\dot{\boldsymbol{X}}_{1}^{\mathrm{T}} \boldsymbol{X}_{1}} \\ {\vdots} \\ {R_{M 1} \dot{R}_{M 1}-\dot{\boldsymbol{X}} V_{M} X_{M}+\dot{\boldsymbol{X}}_{1}^{\mathrm{T}} \boldsymbol{X}_{1}}\end{array}\right] $

$ {\mathit{\boldsymbol{G}}_f} = - 2\left[ {\begin{array}{*{20}{c}} {{{\left( {{{\mathit{\boldsymbol{\dot X}}}_2} - {{\mathit{\boldsymbol{\dot X}}}_1}} \right)}^{\rm{T}}}}&{{{\dot R}_{21}}}&{{{\left( {{\mathit{\boldsymbol{X}}_2} - {\mathit{\boldsymbol{X}}_1}} \right)}^{\rm{T}}}}&{{R_{21}}}\\ {{{\left( {{{\mathit{\boldsymbol{\dot X}}}_3} - {{\mathit{\boldsymbol{\dot X}}}_1}} \right)}^{\rm{T}}}}&{{{\dot R}_{31}}}&{{{\left( {{\mathit{\boldsymbol{X}}_3} - {\mathit{\boldsymbol{X}}_1}} \right)}^{\rm{T}}}}&{{R_{31}}}\\ {}& \vdots &{}&{}\\ {{{\left( {{{\mathit{\boldsymbol{\tilde X}}}_M} - {{\mathit{\boldsymbol{\dot X}}}_1}} \right)}^{\rm{T}}}}&{{{\dot R}_{M{\rm{1}}}}}&{{{\left( {{\mathit{\boldsymbol{X}}_{\rm{M}}} - {\mathit{\boldsymbol{X}}_1}} \right)}^{\rm{T}}}}&{{R_{M{\rm{1}}}}} \end{array}} \right] $

由式(16)、(17)可知时差频差联合估计观测矩阵为：

$ \mathit{\boldsymbol{\varepsilon }} = \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{\varepsilon }}_t}}\\ {{\mathit{\boldsymbol{\varepsilon }}_f}} \end{array}} \right] = \mathit{\boldsymbol{h}} - \mathit{\boldsymbol{G\theta }} $

式中：h=[h_t, h_f]^T，G=[G_t, G_f]^T。

则适应度函数fitness=‖h-Gθ‖。

3 仿真实验

由于无源定位不仅与测量误差有关，还与布站方式有关，在这里只讨论在固定的布站方式下不同测量误差下的仿真，布站方式采用正方形。用标准PSO算法和优化PSO粒子群算法对定位结果的分析。

仿真参数为粒子数目N为200；学习因子c₁、c₂都为1；杂交概率取0.9；杂交池大小比例为0.4；迭代次数为300。站址参数如表 1、2。

实验1  在测时、测频误差的总噪声误差为-20 dB时，标准PSO算法定位结果收敛曲线如图 2。

实验2  在参数相同的情况下，对优化PSO算法的定位结果收敛曲线如图 3。

实验3  在参数与实验1相同的情况下，对标准PSO算法和优化PSO算法收敛速度进行对比，即满足设定好的门限值的迭代次数对比，该门限值为解出位置的相对误差小于1%，并做1 000次蒙特卡洛实验。

实验4  测量误差为-10 dB时多站TDOA/FDOA联合定位的GDOP分析如图 4。

实验5  测量误差为0 dB时多站TDOA/FDOA联合定位的GDOP分析如图 5。

由实验1和实验2可知优化的PSO算法对时差频差联合定位的收敛更快，用的迭代次数少，当迭代次数到达50次时优化的PSO算法得到结果已经完全接近于真实的辐射源位置，优化PSO算法中前期的收敛速度要快于标准PSO算法。实验3证实了在参数相同且种群粒子数目也相同的情况下，优化PSO的TDOA/FDOA联合定位算法在估计出目标位置所需迭代次数短和时间少，由此可看出优化PSO算法的收敛速度要优于标准PSO算法，其中实验3的表 1标准PSO算法迭代次数是优化PSO算法的2倍，一方面由于标准PSO算法出现了局部最优的情况没有达到设定好的门限值而完成所有的迭代次数，另一方面由于标准PSO算法中的惯性权重系数是固定的，这样就使得在进行速度更新时不能达到最优化，从而增加了迭代次数。实验4和实验5分析了TDOA/FDOA联合定位的误差分析，从图 4和图 5可以看出，联合定位的定位精度随着测时、测频误差的增大而变差，定位误差增大，图中等高线相对密集的地方定位效果不好。

4 结论

传统的定位算法一般都是先通过时差方程求出辐射源位置信息，然后再通过FDOA信息对辐射源速度信息进行求解，而本文算法能同时对辐射源位置和速度同时求得。由以上仿真可知优化PSO算法对无源时差频差联合定位收敛速度比传统PSO算法更快，而且得出最优结果的迭代次数更少。由误差分析可知定位精度由测量误差的增加而增大，而且正方形布站存在两个方向的定位盲区。在后续的研究中考虑观测站自定位误差对辐射源定位的影响，并且对算法做进一步的改进来提高定位精度。